Free Access
Issue
A&A
Volume 616, August 2018
Article Number A176
Number of page(s) 9
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201832747
Published online 07 September 2018

© ESO 2018

1. Introduction

We have set up an automated system called NEORANGER that is built on OpenOrb1 (Granvik et al. 2009) and computes asteroid-Earth impact probabilities for new or updated objects on the Near-Earth-Object Confirmation Page (NEOCP2) provided by the Minor Planet Center (MPC). The orbit-computation methods used by NEORANGER are optimized for cases where the amount of astrometry is scarce or it spans a relatively short time span.

The main societal benefit of monitoring for imminent impactors is that it allows for a warning to be sent out in case of an imminent impact that can no longer be prevented. Two well-known examples of damage-causing impacts are the Tunguska (Chyba et al. 1993) and Chelyabinsk (Popova et al. 2013) airbursts. Neither of these two asteroids hit the surface of the Earth intact but instead disrupted in the atmosphere with the resulting shock waves being responsible for the damage on the ground. While the Tunguska event in 1908 destroyed a wide forested area, it did not cause any injuries because it happened in a sparsely populated area. The Chelyabinsk event in 2013, on the other hand, took place in an urban area injuring many residents and causing damage to buildings. The Chelyabinsk asteroid was undetected before its entry into the atmosphere, because it arrived from the general direction of the Sun.

The primary scientific benefit of discovering asteroids prior to their impact with the Earth is that it allows for the characterization of the impactor in space prior to the impact event. Cross-correlating spectroscopic information and the detailed mineralogical information obtained from the meteorites found at the impact location allows us to extend the detailed mineralogical information to all asteroids by using the spectroscopic classification of asteroids as a proxy. Only two asteroids have been discovered prior to entry into the Earth’s atmosphere, 2008 TC3 (Jenniskens et al. 2009) and 2014 AA (Farnocchia et al. 2016a). Both asteroids were small (no more than five meters across) and hence discovered only about 20 h before impact. Whereas hundreds of observations were obtained of 2008 TC3, only seven observations were obtained of 2014 AA. The impact location and time were accurately predicted in advance for 2008 TC3 whereas for 2014 AA the accurate impact trajectory was reconstructed afterwards using both astrometric observations and infrasound data. Meteorite fragments were recovered across the impact location for 2008 TC3 whereas no meteorite fragments could be collected for 2014 AA, which fell into the Atlantic Ocean. Spectroscopic observations in visible wavelengths were obtained for 2008 TC3 suggesting that the object is a rare F-type. Mineralogical analysis of one of the fifteen recovered meteorites showed it to be an anomalous polymict ureilite and thereby established a link between F-type asteroids and ureilites (Jenniskens et al. 2009).

Let us make the simplifying assumption that meteorites will be produced in the impacts of all asteroids larger than one meter in diameter. Brown et al. (2002, 2013) estimate that five-meterdiameter or larger asteroids (that is, similar to 2008 TC3) impact the Earth once every year. Since about 70% of the impacts happen above oceans and roughly 50% of the impactors will come from the general direction of the Sun, our crude estimate is that an event similar to the impact of 2008 TC3 (pre-impact spectroscopic observations and discovery of meteorites) will happen approximately twice in a decade. Since not all of the incoming asteroids are detected, our estimate is to be considered optimistic. This is also supported by the fact that there has not been another 2008 TC3-like event in the past decade. These rates imply that it will take at least a millennium to get a statistically meaningful sample of, say, ten events for each of the approximately 20 Bus-DeMeo spectral classes. However, the rate of events can be increased if the surveys could detect smaller objects prior to their impact.

The Asteroid Terrestrial-impact Last Alert System (ATLAS, Tonry et al. 2018) seeks to provide a warning time in the order of tens of hours for a Chelyabinsk-size object. ATLAS currently consists of two telescopes in the Hawaiian Islands with plans to build additional observatories in other parts of the world. The cameras of these two telescopes have extremely wide fields of view compared to other contemporary near-Earth-object (NEO) survey systems, and they scan the sky in less depth but more quickly producing up to 25 000 asteroid detections per night. The European Space Agency’s (ESA) fly-eye telescope is a similar concept with very large field of view. First light is expected in 2018 and it will be the first telescope in a proposed future European network that would scan the sky for NEOs. We stress that both ATLAS and ESA’s fly-eye telescopes have been funded primarily to reduce the societal risks associated to small asteroids impacting the Earth.

Pushing the detection threshold down to one-meter-diameter asteroids will increase the discovery rate for Earth impactors by a factor of 40. A factor of five in physical size corresponds to three magnitudes and such an improvement will, optimistically, be provided by the Large Synoptic Survey Telescope (LSST) with first light in 2023. Since events similar to 2008 TC3 will still be rare and the timescales short, it is of the utmost importance to identify virtually all impactors and do that as early as possible, that is, immediately after their discovery, to allow for spectroscopic follow-up. Another challenge will be to obtain useful spectroscopic data of small, fast-moving impactors; obtaining visible spectroscopy of 2008 TC3 was a borderline case with a four-meter-class telescope.

An alternative approach to determine the compositional structure of the asteroid belt is to focus on meteorite falls and use their predicted trajectories prior to the impact to identify the most likely source region in the asteroid belt. Granvik & Brown (in prep.) use the observed trajectories of 25 meteorite-producing fireballs published to date to associate meteorite types to their most likely escape routes from the asteroid belt and the cometary region. We note that while this approach is indirect, because it lacks the spectroscopic observations prior to the impact, it results in substantially higher statistics per unit time.

Long-term (timescales ranging from months to centuries) asteroid impact monitoring is carried out by Jet Propulsion Laboratory’s (JPL’s) Sentry3 system (Chamberlin et al. 2001) and the University of Pisa’s close-approach monitoring system CLOMON24 (Chesley & Milani 1999). These automated systems are based on a similar algorithm (Milani et al. 2002, 2005; Farnocchia et al. 2015b) and run in parallel. The results are continuously compared to identify problematic cases that have to be scrutinized in greater detail by human operators.

Whereas the long-term impact monitoring is well-established, dedicated short-term (from hours to weeks) impact monitoring has only recently started receiving serious attention. Farnocchia et al. (2015a) implemented an automated system (Farnocchia et al. 2016b) to regularly analyze the objects on the MPC’s NEOCP with systematic ranging. Spoto et al. (2018) have developed a short-arc orbit computation method that also analyzes the objects on the NEOCP with systematic ranging. Also NEORANGER relies on the observations on the NEOCP, but the implementation and the orbit-computation method are different. As in the case for long-term impact monitoring, comparing the results of two or more systems doing the same analysis but using different techniques provides assurance that the results are valid. The very time-critical science described above underscores the need for accurate and complete monitoring. In addition to scientific utility, small impactors such as 2008 TC3 and 2014 AA are valuable in comparing the predicted and actual impact location and time. These small, non-hazardous impactors thus allow the community to verify that the impact-monitoring systems work correctly.

NEORANGER also monitors for other types of interesting objects such as asteroids temporarily captured by the Earth-Moon system. Granvik et al. (2012) predict that there is a one-meter-diameter or larger NEO temporarily orbiting the Earth at any given time. The population of temporarily-captured natural Earth satellites (NES) are challenging to discover with current surveys because they are small in size and move fast when at a sufficiently small distance (Bolin et al. 2014). Asteroid 2006 RH120 (Kwiatkowski et al. 2009) is to date the only NES confirmed not to be a man-made object. It is a few meters in diameter and therefore less challenging to discover compared to smaller objects.

In what follows we will first describe the real data available through NEOCP as well as the simulated data used for testing NEORANGER. Then we describe the NEORANGER system and the overall characteristics of the numerical methods used (with detailed descriptions to be found in the appendix). Finally, we present a selection of results and discuss their implications, and end with some conclusions and future avenues for improvement.

2. Data

The MPC is a clearinghouse for asteroid and comet astrometry obtained by observers worldwide. Based on an automatically computed probability of a discovered moving object being potentially a new NEO, it is placed on MPC’s NEOCP for confirmation through follow-up observations. Systems such as NEORANGER and JPL’s Scout system5 monitor the NEOCP for astrometry on new discoveries.

Once the orbital solution for a new object has been sufficiently well constrained to secure its orbital classification, the MPC issues a Minor Planet Electronic Circular (MPEC) where the discovery is published, and the object is removed from the NEOCP. Following the release of an MPEC long-term (typically up to 100 yr into the future) impact-probability computations are undertaken by Sentry and CLOMON2.

2.1. Objects on the NEOCP

As the primary test data set, we use 2339 observation sets of 695 objects on the NEOCP between August 19 and October 12, 2015 (eight weeks). We use the term “observation set” to refer to astrometric observations for a single object posted on the NEOCP within the last half an hour, that is, NEORANGER checks the NEOCP for updates every half an hour by default.

Of the 695 objects that appeared on the NEOCP, 260 objects eventually received an MPEC. The remaining 435 objects were not NEOs or other objects on peculiar orbits, did not receive enough follow-up observations to reveal their true nature, were linked to previously discovered asteroids, or were found to be artificial satellites or image artifacts. Zero percent (36%) of the 260 objects that received an MPEC (of the 435 objects that did not receive an MPEC) have only one observation set. Zero percent (36%) have two observation sets, that is, the object is updated once with more observations. Two percent (16%) have three observation sets and 51% (37%) four observation sets. A large part of the observations in our test data set were obtained by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) 1 (F51) as well as the Catalina Sky Survey (703) and the Mt. Lemmon Survey (G96). The first observation set for an object contains observations from only one observatory. Seventy-seven percent of the 695 objects have at least two observation sets, in which case for 91% of the objects the first two observation sets come from different observatories. Fifty-five percent of the 695 objects have at least three observation sets, in which case for 91% of the objects the first three observation sets come from three observatories. An observation set typically contains three or four observations. The time span of an observation set for Pan-STARRS 1 is typically 30–60 min, and for the Catalina Sky Survey and the Mt. Lemmon Survey 15–30 min. The time between two observation sets ranges from about 30 min to two days. The total time span when including up to four observation sets is typically between 2 h and three days.

2.2. Simulated Earth impactors

Apart from the two known Earth impactors mentioned in the introduction, there are no real impactors for testing neoranger. Therefore we tested NEORANGER with 50 simulated objects that were randomly picked from the simulated sample generated by Vereš et al. (2009). We generated synthetic astrometry by using the location and the typical, but simplified, Saturday, June 9, 2018 at 11:02 am cadence of Pan-STARRS 1. For 30 simulated impactors we generated three astrometric observations with a 0.01 day interval between the observations about five days prior to the impact with the Earth. Then we generated three more astrometric observations for the next night, that is, four days before the impact. For the remaining 20 simulated impactors we did the same but generated four astrometric observations instead of three.

2.3. Geocentric objects

Similarly to Earth impactors, we have very few known NESs. Therefore we test the system’s capability to identify NESs with one asteroid temporarily captured by the Earth, 2006 RH120, and also the space observatory Spektr-R. 2006 RH120 was a natural temporarily-captured orbiter, which orbited the Earth from July 2006 until July 2007. For 2006 RH120 we used 17 observations between September 16 and 17, 2006 spanning 24 h divided into four tracklets. The heliocentric elements at the end of the epoch used are a = 1.04 AU, e = 0.031, and i = 1.43° and the geocentric elements are a = 0.013 AU, e = 0.56, and i = 64.2°. Spektr-R (or RadioAstron) is a Russian Earth-orbiting space observatory (Kardashev et al. 2013). For Spektr-R we used eight observations on April 4, 2016 spanning 4 h divided into three tracklets. The geocentric elements are a = 180 000 km, e = 0.93, and i = 38.8° at the end of the epoch used.

3. Methods

3.1. Statistical ranging

In statistical orbital ranging (Virtanen et al. 2001; Muinonen et al. 2001), thetypicallynon-Gaussianorbital-elementprobabilitydensity is examined using Monte Carlo selection of orbits in orbital-element space using the Bayesian formalism. In practice, the statistical ranging method proceeds in the following way. Two observations are selected from the data set and random deviates are added to the four plane-of-sky coordinates to mimic observational noise. Next, topocentric distances are generated for the two observation dates. The noisy sky-plane coordinates and the two topocentric distances are transformed, by accounting for the observatory’s heliocentric coordinates at the observation dates, into two heliocentric position vectors that are mapped into the phase space of a Cartesian state vector, which is equivalent to a set of orbital elements (e.g., Keplerian elements). The proposed orbit is then used to compute ephemerides for all observation dates. In Monte Carlo (MC) ranging the proposed orbit is accepted if it produces acceptable sky-plane residuals and the Δχ2 is small enough with respect to the until-then best-fit orbit. The two other ranging variants using Markov chains require a burn-in phase, which ends when, for random walk (Muinonen et al. 2016), the first orbit with acceptable Δχ2 is found, and for Markov-Chain Monte Carlo (MCMC, Oszkiewicz et al. 2009), the first orbit with acceptable residuals at all dates is found. We use a constant prior probability-density function (PDF) with Cartesian state vectors.

The different variants (MC, MCMC, and random-walk) of the statistical orbital ranging method (Appendix A) used here have been implemented in the open-source orbit-computation software package OpenOrb. The ranging method in general has been tested with several large-scale applications, demonstrating its wide applicability to various observational data (see Virtanen et al. 2015, and references therein).

The systematic ranging used by Scout (Farnocchia et al. 2015a) scans a grid in the space of topocentric range ρ and range rate ρ̇. For each grid point a fit in right ascension and declination as well as their time derivatives (α, δ, α̇, δ̇) is computed by minimizing a cost function created from the observed-computed residuals and, combined with the topocentric range and range rate, these six are equivalent to a complete set of orbital elements. Farnocchia et al. (2015a) also tested different priors to be used with their analysis and concluded that a uniform prior produces the best results. The choice of an adequate prior distribution is critical to ensuring that impacting asteroids are identified as early as possible. In what follows we will use the uniform prior because that has been our default for many years and now also Farnocchia et al. (2015a) have shown that it is indeed a good choice. An alternative implementation of the systematic ranging was recently proposed by Spoto et al. (2018).

3.2. NEORANGER system description

The NEORANGER system evaluates the Earth-impact probability and various other characteristics of objects on the NEOCP. The automated system works in general as follows. NEORANGER checks the NEOCP every 30 min for new objects and old objects with new data. We note that the frequency of checks is a tunable parameter and can be changed, but we have empirically found out that 30 min is a reasonable update frequency for our current monitoring purposes.

For each observation set the different ranging methods are run (with a timeout of 20 min) in the following order until the requested number of sample orbits are computed: Monte Carlo (MC), Markov-Chain Monte Carlo (MCMC), and random-walk ranging (Appendix A). Previous orbit solutions are provided as input to ranging to constrain the initial range bounds and speed up the computations. For each observation set that results in the requested number of accepted sample orbits, the orbital-element PDF is propagated forward in time for seven days to derive the impact parameter distribution. We note that also the length of the forward propagation is a configuration parameter that can be changed if deemed necessary. The computing time required increases roughly linearly with the length of the forward propagation. The time needed to identify impactors can be reduced by dividing the integrations according to a logarithmic sequence, say, one, three, ten, and 30 days, so that, for example, the alert for an asteroid impacting within a day can be sent out after the first integration rather than waiting for the 30 day integration to finish. The impact probability is computed as the weighted fraction of impact orbits. Similarly, the probability of the object being geocentric is computed by transforming the heliocentric orbital elements to geocentric orbital elements, and computing the ratio of those orbits that have geocentric e < 1.

The NEORANGER system sends out a notification by e-mail of all cases where the impact probability or the probability for an object being captured in an orbit about the Earth is greater than 10−4. We will next explain how this critical value and the required number of sample orbits have been determined.

4. Results

4.1. Configuration parameters

We used conservative values for the many adjustable configuration parameters in OpenOrb to be less sensitive to outliers in the first data sets and to allow for the widest possible orbital-element distribution. For the astrometric uncertainty we therefore assumed one arcsecond for all observatories. The configuration parameter that has the largest influence on the computing time is the required number of sample orbits. Since the number of sample orbits also directly affects the reliability of the results, we decided to carry out tests to find a suitable number. After 50 000 orbits the probabilities start to converge when gradually increasing the number of required sample orbits towards 100 000 (Fig. 1). We therefore fixed the number of required sample orbits to 50 000, because this is a good compromise between computing time and the consistency of the results.

thumbnail Fig. 1.

Impact probabilities for 398 objects with 100 000 and 50 000 orbits as the number of required sample orbits.

Open with DEXTER

4.2. Earth impactors

We tested the system with the two known Earth impactors (2008 TC3 and 2014 AA) and 50 objects from the simulated sample generated by Vereš et al. (2009). For the latter we generated synthetic astrometry for two consecutive nights as explained in Sect. 2.2. To compute the sample orbits the MC ranging was run twice: first with the observations of the first night, and then with the observations of both nights initializing the analysis with the orbits from the first ranging. The orbital-element probability density function resulting from the ranging is propagated forward in time as explained in Sect. 3.2. We found a strong correlation between the impact probability and the number of impacting orbits for the simulated impactors (Fig. 2). This implies that high impact probabilities are based on a statistically meaningful number of sample orbits, and are thus unlikely to be false alarms. We note that we do not here account for the fact that the typical false alarm is caused by a so-called outlier observation, an astrometric observation with the corresponding residual substantially larger than the expected astrometric uncertainty. The smallest impact probability is 10%, ten objects have a probability less than 50%, and 32 objects have a probability greater than 90%. For the nights with four astrometric observations the impact probabilities do not increase systematically compared to the nights with three observations.

thumbnail Fig. 2.

For 50 simulated impactors, the impact probability versus the fraction of impacting orbits. The solid black line is included for reference and corresponds to a 100% correlation between the fraction of impacting orbits and the impact probability.

Open with DEXTER

When testing NEORANGER with the two known impactors, 2008 TC3 and 2014 AA, we use two tracklets built from the first data sets available. For 2008 TC3 we use the first observations obtained by R. Kowalski: for the first tracklet four observations and for the second tracklet three additional observations. For 2014 AA we use three observations for the first tracklet and the remaining four observations available for the second tracklet. For the first tracklet the ranging takes 1 min and for the second tracklet 10 min for both objects. With the addition of the second tracklet the propagation takes a considerably longer time, an hour for 2008 TC3 and an hour and a half for 2014 AA), because due to many impacting orbits the integration steps become smaller. For both objects the impact probabilities are negligible (less than one in ten billion) when using only the first tracklet with a uniform prior, but increase to 81%–60%, respectively, with Jeffreys’ prior. The apparent discrepancy with the results obtained with a uniform prior by Farnocchia et al. (2015a) will be studied in detail in a forthcoming paper. With the addition of the second tracklet the impact probabilities increase to 92% for 2008 TC3 and to 80% for 2014 AA. The two-dimensional (2D) marginal probability densities for semi-major axis and eccentricity for 2008 TC3 and 2014 AA are shown in Fig. 3. The PDFs resulting from the ranging method correctly cover the region where the “true” orbit based on all available data (marked with a star) resides.

thumbnail Fig. 3.

Two-dimensional marginal probability densities for semi-major axis and eccentricity for the impactors 2008 TC3 (epoch 54745.8 MJD) and 2014 AA (epoch 56658.3 MJD) based on the first seven observations for both objects. The “true” orbital elements for the chosen epoch based on all available data as reported by JPL’s HORIZONS system (https://ssd.jpl.nasa.gov/horizons.cgi) are marked with a star. We note that the estimate for 2014 AA also includes infrasound data (Farnocchia et al. 2016a).

Open with DEXTER

In Fig. 4 the red histogram on the right represents 2008 TC3, 2014 AA, and 50 simulated objects using two observation sets from consecutive nights. The green histogram represents 50 simulated objects using only one observation set. Hence it is clear that the substantially increased impact probability when adding a single data set to the discovery data set is typical for all impactors rather than being a random occurence for 2008 TC3 and 2014 AA.

thumbnail Fig. 4.

Impact probabilities for NEOCP objects during eight weeks of continuous monitoring (blue), for 50 simulated impactors using one observation set (green), and for 2008 TC3, 2014 AA, and 50 simulated impactors using two observation sets from consecutive nights (red).

Open with DEXTER

4.3. Geocentric objects

In addition to potential Earth-impacting objects our system also monitors for Earth’s natural temporarily-captured satellites, that is, objects on elliptic, geocentric orbits (hereafter just geocentric orbits). Here we present results of tests with 2006 RH120, an asteroid temporarily captured by the Earth, and the space observatory Spektr-R.

For 2006 RH120 the first tracklet (four observations during 3 min) produces only a negligible probability for the object being geocentric, the second tracklet (nine observations during 7 h) a 23% probability, the third tracklet (14 observations during 7 h) a 97% probability, and the fourth tracklet (17 observations during 24 h) a 100% probability. For Spektr-R the first tracklet (four observations during 1 h) gives a 24% probability, the second tracklet (six observations during 3 h) a 97% probability, and the third tracklet (eight observations during 4 h) a 100% probability. From these two examples we see that our system recognizes a geocentric object only if both the time span and number of observations are sufficient.

4.4. NEOCP objects

Having verified that NEORANGER can identify Earth-impacting asteroids as well as objects on elliptical, geocentric orbits, we now present results of eight weeks of continuous operations between August 18 and October 12 2015 for objects on the NEOCP. For each object NEORANGER calculates the probability that the object will impact the Earth within a week as well as the probability that the object is on an elliptic, geocentric orbit. In addition, we compute an estimate for the perihelion distance and NEO class. We compare these results with those given by MPC.

We use 2339 observation sets of 695 objects on the NEOCP. To derive the impact parameter distribution the ranging is run for each observation set and the resulting orbital element probability density function is propagated forward in time for seven days as explained in Sect. 3.2. The ranging for the first observation set takes on average 1 min with root-mean-square (RMS) deviation of 3 min, for the second observation set 3 min (RMS 5 min), for the third observation set 6 min (RMS 7 min), and for the fourth observation set 7 min (RMS 8 min). The ranging gradually slows down or fails for observation sets with larger time spans and number of observations. If the ranging fails it is mostly because the timeout of 20 min is exceeded or because no acceptable sample orbits are found (e.g., because the residuals are too large compared to the assumed astrometric uncertainty). The time required for the propagation stays typically around 1.5 min.

A non-zero impact probability is computed for 424 observation sets (and in accordance with Virtanen & Muinonen (2006) almost always for the first observation set of the object). The blue histogram in Fig. 4 shows the distribution of all non-zero impact probabilities computed by NEORANGER during the eight weeks of continuous monitoring. The ten highest probabilities in the blue histogram fall between 7.7 × 10−10 and 2.3 × 10−8. There is a clear gap between, on one hand, the group formed by the 50 simulated impactors and the two real impactors and, on the other hand, the group formed by the NEOCP objects and the simulated impactors with only one observation set. It thus seems that impactors are clearly identified already based on only two observation sets, and NEORANGER produces very few, if any, false alarms.

A non-zero probability for the object being on an elliptic, geocentric orbit is computed for 727 observation sets. Some fraction of these objects correspond to artificial Earth-orbiting satellites. In particular, NEORANGER has computed a non-zero probability for a geocentric orbit for all artificial objects that have appeared on NEOCP and that have been identified as such on the various discussion forums.

The evolution of the perihelion-distance estimate with the increasing number of observation sets for the 260 objects that received an MPEC is shown in Fig. 5. The error bars define the 1σ confidence interval. The accuracy increases with each new observation set as expected. Figure 6 shows examples of the perihelion-distance distributions for two objects. We note the substantially non-Gaussian distribution in the top left plot compared to the Gaussian distribution in the bottom left plot.

thumbnail Fig. 5.

Perihelion distance as estimated by NEORANGER and by MPC. The MPC estimate is based on all available data whereas the neoranger estimate only includes parts of it as described below each plot. The solid red line is included for reference and corresponds to a 100% correlation between the perihelion distance as estimated by NEORANGER and by MPC.

Open with DEXTER

thumbnail Fig. 6.

Two example distributions of the perihelion distance. The red vertical dashed line corresponds to the maximum likelihood and the blue vertical lines define the 1σ confidence interval.

Open with DEXTER

NEORANGER predicts the NEO class (Amor, Apollo, Aten, Aethra/Mars-crosser) correctly (i.e., NEORANGER estimates the class to be the one given by MPC with over 50% confidence) for 40% of the 260 objects in our test data set when using only the first observation set and 63% when using two observation sets (Fig. 7). In the rightmost bar in the plot for two tracklets we see that the correct class is predicted for 0.377 × 260 = 98 objects with at least a 90% probability.

thumbnail Fig. 7.

Fractional distribution of the probabilities (as estimated by NEORANGER) to belong to the class reported by the MPC when using one or two observation sets.

Open with DEXTER

In rare cases (one object out of a thousand based on 13 months of operations) NEORANGER erroneously computes a high impact probability based on a small number of impacting orbits, as in the MCMC ranging run for 2017 JK2 (Fig. 8). In that case, only 12 impacting orbits produce an impact probability of 2% although all the impacting orbits have large χ2 values. These impacting orbits also correspond to very small topocentric distances at the orbit computation epoch, which is typically the mid-point of the observational time span. These impacting orbits thus comprise the low-likelihood tail of the orbit distribution. In MCMC ranging the weight of a sample orbit is computed based on the repetitions, that is, the inability of the proposed orbit to be accepted when starting from the sample orbit in question. A viable scenario therefore is that, in some rare cases, the sample orbits in the tail get an unrealistically large number of repetitions, because there is essentially just one narrow path towards higher-likelihood orbits and the algorithm requires an unusually long time to find it. Since a large fraction of these close-in orbit solutions typically impact the Earth, the impact probability is also substantial despite the fact that the solutions have low likelihoods based on their χ2 values. The other two examples in Fig. 8 are the real impactor 2008 TC3 and the simulated impactor Si000003 from Vereš et al. (2009). For 2008 TC3 we get an impact probability of 98% with 48 293 impacting orbits (out of 50 000 orbits), and for Si000003 an impact probability of 24% with 11 691 impacting orbits.

thumbnail Fig. 8.

Distribution of PDFs and χ2 values for 2008 TC3, Si000003, and 2017 JK2 resulting from one MCMC ranging run. The red symbols represent the impacting orbits.

Open with DEXTER

5. Conclusion

While NEORANGER correctly predicts the impacts of the two asteroids discovered before impacting the Earth, 2008 TC3 and 2014 AA, the impact probabilities for objects on the NEOCP computed by NEORANGER are always less than one in ten million, as it should be for non-impactors typically found on NEOCP. In addition, NEORANGER rarely produces false alarms of impacting asteroids. As we gain experience from real objects we will be able to apply a threshold for an impact probability that should trigger follow-up efforts or, at least, close scrutiny of the data.

We note that, for the two known impactors (2008 TC3 and 2014 AA) and 50 simulated impactors, a significant impact probability is achieved only using two tracklets while using only one tracklet results in a negligible impact probability. Asaconsequence our system might not be able to identify an impactor based on the very first data set available. It may be possible to increase the impact probability computed for the first observation set by using a uniform prior in polar coordinates (Farnocchia et al. 2015a). We leave the implementation and testing of the alternative prior(s) to future works that should also provide a quantitative comparison of the existing short-term impact-monitoring systems. Such a comparison would provide end-users with a fact-based list of pros and cons of the different systems.


Acknowledgments

We thank Federica Spoto (reviewer) and Davide Farnocchia for comments and suggestions that helped improve the paper. OS acknowledges funding from the Waldemar von Frenckell foundation. MG acknowledges funding from the Academy of Finland (grant #299543).

References

  1. Bolin, B., Jedicke, R., Granvik, M., et al. 2014, Icarus, 241, 280 [NASA ADS] [CrossRef] [Google Scholar]
  2. Brown, P., Spalding, R. E., ReVelle, D. O., Tagliaferri, E., & Worden, S. P. 2002, Nature, 420, 294 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brown, P. G., Assink, J. D., Astiz, L., et al. 2013, Nature, 503, 238 [NASA ADS] [Google Scholar]
  4. Chamberlin, A. B., Chesley, S. R., Chodas, P. W., et al. 2001, in AAS/Division for Planetary Sciences Meeting Abstracts, BAAS, 33, 1116 [Google Scholar]
  5. Chesley, S. R., & Milani, A. 1999, in AAS/Division for Planetary Sciences Meeting Abstracts, 31, 28.06 [Google Scholar]
  6. Chyba, C. F., Thomas, P. J., & Zahnle, K. J. 1993, Nature, 361, 40 [NASA ADS] [CrossRef] [Google Scholar]
  7. Farnocchia, D., Chesley, S. R., & Micheli, M. 2015a, Icarus, 258, 18 [NASA ADS] [CrossRef] [Google Scholar]
  8. Farnocchia, D., Chesley, S. R., Milani, A., Gronchi, G. F., & Chodas, P. W. 2015b, Orbits, Long-Term Predictions, Impact Monitoring, eds. P. Michel, F. E. DeMeo, & W. F. Bottke, 815 [Google Scholar]
  9. Farnocchia, D., Chesley, S. R., Brown, P. G., & Chodas, P. W. 2016a, Icarus, 274, 327 [NASA ADS] [CrossRef] [Google Scholar]
  10. Farnocchia, D., Chesley, S. R., & Chamberlin, A. B. 2016b, in AAS/Division for Planetary Sciences Meeting Abstracts, 48, 305.03 [Google Scholar]
  11. Granvik, M., Virtanen, J., Oszkiewicz, D., & Muinonen, K. 2009, Meteor. Planet. Sci., 44, 1853 [NASA ADS] [CrossRef] [Google Scholar]
  12. Granvik, M., Vaubaillon, J., & Jedicke, R. 2012, Icarus, 218, 262 [NASA ADS] [CrossRef] [Google Scholar]
  13. Jenniskens, P., Shaddad, M. H., Numan, D., et al. 2009, Nature, 458, 485 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  14. Kardashev, N. S., Khartov, V. V., Abramov, V. V., et al. 2013, Astron. Rep., 57, 153 [NASA ADS] [CrossRef] [Google Scholar]
  15. Kwiatkowski, T., Kryszczyńska, A., Polińska, M., et al. 2009, A&A, 495, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Milani, A., Chesley, S. R., Chodas, P. W., & Valsecchi, G. B. 2002, in Asteroid Close Approaches: Analysis and Potential Impact Detection, eds. W. F. Bottke, Jr., A. Cellino, P. Paolicchi, & R. P. Binzel, 55 [Google Scholar]
  17. Milani, A., Chesley, S. R., Sansaturio, M. E., Tommei, G., & Valsecchi, G. B. 2005, Icarus, 173, 362 [NASA ADS] [CrossRef] [Google Scholar]
  18. Muinonen, K., & Bowell, E. 1993, Icarus, 104, 255 [NASA ADS] [CrossRef] [Google Scholar]
  19. Muinonen, K., Virtanen, J., & Bowell, E. 2001, Celest. Mech. Dyn. Astron., 81, 93 [NASA ADS] [CrossRef] [Google Scholar]
  20. Muinonen, K., Fedorets, G., Pentikäinen, H., et al. 2016, Planet. Space Sci., 123, 95 [NASA ADS] [CrossRef] [Google Scholar]
  21. Oszkiewicz, D., Muinonen, K., Virtanen, J., & Granvik, M. 2009, Meteor. Planet. Sci., 44, 1897 [NASA ADS] [CrossRef] [Google Scholar]
  22. Popova, O. P., Jenniskens, P., Emel’yanenko, V., et al. 2013, Science, 342, 1069 [NASA ADS] [CrossRef] [Google Scholar]
  23. Spoto, F., Del Vigna, A., Milani, A., et al. 2018, A&A, 614, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, PASP, 130, 064505 [NASA ADS] [CrossRef] [Google Scholar]
  25. Vereš, P., Jedicke, R., Wainscoat, R., et al. 2009, Icarus, 203, 472 [NASA ADS] [CrossRef] [Google Scholar]
  26. Virtanen, J., & Muinonen, K. 2006, Icarus, 184, 289 [NASA ADS] [CrossRef] [Google Scholar]
  27. Virtanen, J., Muinonen, K., & Bowell, E. 2001, Icarus, 154, 412 [NASA ADS] [CrossRef] [Google Scholar]
  28. Virtanen, J., Fedorets, G., Granvik, M., Oszkiewicz, D. A., & Muinonen, K. 2015, IAU General Assembly, 22, 56438 [NASA ADS] [Google Scholar]

Appendix A: Statistical orbital ranging

A.1. Inverse problem of orbit computation

Within the Bayesian framework, the orbital-element PDF is(A.1)

where the prior PDF ppr is constant and the observational error PDF p is evaluated for the observed–computed sky-plane residuals ΔΨ(P) and is typically assumed to be Gaussian (Muinonen and Bowell 1993). The parameters P describe the orbital elements of an asteroid at a given epoch t0. For Keplerian orbital elements P = (a, e, i, Ω, ω, M0)T. For Cartesian elements, P = (X, Y, Z, Ẋ, Ẏ, Ż)T, where, in a given reference frame, the vectors (X, Y, Z)T and (Ẋ, Ẏ, Ż)T denote the position and velocity, respectively.

A.2. Statistical ranging

In statistical orbital ranging (Virtanen et al. 2001; Muinonen et al. 2001), the orbital-element PDF is examined using Monte Carlo selection of orbits in orbital-element space in the following way:

  • Select two observations (usually the first and the last) from the data set.

  • To vary the topocentric coordinates R.A. and Dec., and the topocentric range of the two observations introduce six uniform random deviates mimicking uncertainties in the coordinates.

  • Compute new Cartesian positions at the dates.

  • For the proposed orbit, compute candidate orbital elements, R.A. and Dec. for all observation dates, and Δχ2.

  • Let(A.2)

    where Λ is the covariance matrix for the observational errors and P 0 specifies a reference orbital solution. Notice that, for linear models and Gaussian PDFs, the definition of Eq. (A.2) yields the well-known result(A.3)

    where P0 denotes the linear least-squares orbital solution. In MC ranging proposed orbits are accepted if they produce acceptable sky-plane residuals: the Δχ2 value of the residuals is below a given threshold (we use 30), and if the residuals at all observation dates are smaller than given cutoff values (for Δαmax and Δδmax we use 4 arcsec). For the acceptance criteria of the other ranging variants see Sects. A.3 and A.4.

  • If accepted, assign a statistical weight based on χ2(P) describing the orbital-element PDF of the sample orbit:(A.4)

  • The procedure is repeated up to 107 times resulting in up to 50 000 accepted sample orbits. Both of these are adjustable parameters and their current values have been found empirically (see Sect. 4.1).

A.3. MCMC ranging

Markov-chain Monte Carlo (MCMC) ranging (Oszkiewicz et al. 2009) is initiated with the selection of two observations from the full set of observations: typically, the first and the last observation are selected, denoted by A and B. Orbital-element sampling is then carried out with the help of the corresponding topocentric ranges (ρA, ρB), R.A.s (αA, αB), and Dec.s (δA, δB). These two spherical positions, by accounting for the light time, give the Cartesian positions of the object at two ephemeris dates. The two Cartesian positions correspond to a single, unambiguous orbit passing through the positions at the given dates.

For sampling the orbital-element PDF we utilize the Metropolis-Hastings (MH) algorithm. The MH algorithm is based on the computation of the ratio ar:(A.5)

where Pj and Qj denote the current orbital elements and spherical positions, respectively, in a Markov chain, and the primed symbols, P and Q, denote their proposals. The proposal PDFs from Pj to P and from Qj to Q (t stands for transition) are, respectively, pt(P; Pj) and pt(Q; Qj).

The proposal P is transformed to the space of two topocentric spherical positions resulting in a multivariate Gaussian proposal PDF pt(Q; Qj). This transformation introduces Jacobians Jj and J′:(A.6)

The ratio ar in Eq. (A.5) simplifies into its final form because the proposal PDFs pt(Qj; Q) and pt(Q; Qj) are symmetric.

The proposed elements P are accepted with the probability of min(1, ar):(A.7)

After a number of transitions in the so-called burn-in phase, the Markov chain, in the case of success, converges to sample the target PDF pp when the first orbit with acceptable residuals at all dates is found.

The posterior distribution is proportional to the number of repetitions of a given orbit r(P) rather than on the χ2(P) values:(A.8)

That is, the difficulty to find an acceptable P increases the probability of the current orbit P.

A.4. Random walk ranging

Instead of MCMC ranging, it can be advantageous to sample in the entire phase-space regime below a given χ2(P) level, assigning weights on the basis of the a posteriori probability density value and the Jacobians presented above (Muinonen et al. 2016; cf., Virtanen et al. 2001; Muinonen et al. 2001).

MCMCranging canbemodifiedfor whatwecall random-walk ranging of the phase space within a given Δχ2 level. First, we assign a constant, nonzero PDF for the regime of acceptable orbital elements and assign a zero or infinitesimal PDF value outside the regime. MCMC sampling then returns a set of points that, after convergence to sampling the phase space of acceptable orbital elements (i.e., the burn-in phase ends when the first orbit with acceptable Δχ2 is found), uniformly characterizes the acceptable regime. Second, we assign the posterior PDF values pp as the weights wj for the sample orbital elements Pj. In the case of ranging using the topocentric spherical coordinates, the weightswj need to be further divided by the proper Jacobian value Jj.

In random-walk ranging, uniformly sampling the phase space of the acceptable orbital elements, the final weight factor for the sample elements Pj in the Markov chain is(A.9)

We note that the same orbit can repeat itself in the chain, analogously to MCMC ranging.

All Figures

thumbnail Fig. 1.

Impact probabilities for 398 objects with 100 000 and 50 000 orbits as the number of required sample orbits.

Open with DEXTER
In the text
thumbnail Fig. 2.

For 50 simulated impactors, the impact probability versus the fraction of impacting orbits. The solid black line is included for reference and corresponds to a 100% correlation between the fraction of impacting orbits and the impact probability.

Open with DEXTER
In the text
thumbnail Fig. 3.

Two-dimensional marginal probability densities for semi-major axis and eccentricity for the impactors 2008 TC3 (epoch 54745.8 MJD) and 2014 AA (epoch 56658.3 MJD) based on the first seven observations for both objects. The “true” orbital elements for the chosen epoch based on all available data as reported by JPL’s HORIZONS system (https://ssd.jpl.nasa.gov/horizons.cgi) are marked with a star. We note that the estimate for 2014 AA also includes infrasound data (Farnocchia et al. 2016a).

Open with DEXTER
In the text
thumbnail Fig. 4.

Impact probabilities for NEOCP objects during eight weeks of continuous monitoring (blue), for 50 simulated impactors using one observation set (green), and for 2008 TC3, 2014 AA, and 50 simulated impactors using two observation sets from consecutive nights (red).

Open with DEXTER
In the text
thumbnail Fig. 5.

Perihelion distance as estimated by NEORANGER and by MPC. The MPC estimate is based on all available data whereas the neoranger estimate only includes parts of it as described below each plot. The solid red line is included for reference and corresponds to a 100% correlation between the perihelion distance as estimated by NEORANGER and by MPC.

Open with DEXTER
In the text
thumbnail Fig. 6.

Two example distributions of the perihelion distance. The red vertical dashed line corresponds to the maximum likelihood and the blue vertical lines define the 1σ confidence interval.

Open with DEXTER
In the text
thumbnail Fig. 7.

Fractional distribution of the probabilities (as estimated by NEORANGER) to belong to the class reported by the MPC when using one or two observation sets.

Open with DEXTER
In the text
thumbnail Fig. 8.

Distribution of PDFs and χ2 values for 2008 TC3, Si000003, and 2017 JK2 resulting from one MCMC ranging run. The red symbols represent the impacting orbits.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.