NEOPROP: a NEO Propagator for Space Situational Awareness
V. Zuccarelli, D. Bancelin, S. Weikert, W. Thuillot, D. Hestroffer, C.Y. Valles, D. Koschny
NNEOPROP: A NEO PROPAGATOR FOR SPACE SITUATIONAL AWARENESS FOR THE 6 TH IAASS CONFERENCE
Valentino Zuccarelli (1) , David Bancelin (2) ,Sven Weikert (3) , William Thuillot (4) , Daniel Hestroffer (5) , Celia YabarValles (6) , Detlef Koschny (7) (1)
Astos Solutions GmbH, Grund 1, 78089 Unterkirnach (Germany), Email: [email protected] (2)
IMCCE-Observatoire de Paris, Av. Denfert-Rochereau 77, 75014 Paris (France), Email: [email protected] (3)
Astos Solutions GmbH, Grund 1, 78089 Unterkirnach (Germany), Email: [email protected] (4)
IMCCE-Observatoire de Paris, Av. Denfert-Rochereau 77, 75014 Paris (France), Email: [email protected] (5)
IMCCE-Observatoire de Paris, Av. Denfert-Rochereau 77, 75014 Paris (France), Email: [email protected] (6)
ESA-ESTEC, Keplerlaan 1, 2201 Noordwijk, (The Netherlands), Email: [email protected] (7)
ESA-ESTEC, Keplerlaan 1, 2201 Noordwijk, (The Netherlands), Email: [email protected]
ABSTRACT
The overall aim of the Space Situational Awareness(SSA) Preparatory Programme is to support theEuropean independent utilisation of and access to spacefor research or services, through providing timely andquality data, information, services and knowledgeregarding the environment, the threats and thesustainable exploitation of the outer space surroundingour planet Earth. The SSA system will comprise threemain segments: Space Weather (SWE) monitoring and forecast Near-Earth Objects (NEO) survey and follow-up Space Surveillance and Tracking (SST) of man-made space objectsCurrently, there are over 600.000 asteroids known in ourSolar System, where more than 9.500 of these areNEOs. These could potentially hit our planet anddepending on their size could produce considerabledamage. For this reason NEOs deserve active detectionand tracking efforts.The role of the SSA programme is to provide warningservices against potential asteroid impact hazards,including discovery, identification, orbit prediction andcivil alert capabilities. ESA is now working to develop aNEO Coordination Centre which will later evolve into aSSA-NEO Small Bodies Data Centre (SBDC), locatedat ESA/ESRIN, Italy. The Software prototype developedin the frame of this activity may be later implemented asa part of the SSA-NEO programme simulators aimed atassessing the trajectory of asteroids.There already exist different algorithms to predict orbitsfor NEOs. The objective of this activity is to come upwith a different trajectory prediction algorithm, whichallows an independent validation of the currentalgorithms within the SSA-NEO segment (e.g.NEODyS, JPL Sentry System). The key objective of this activity was to design,develop, test, verify, and validate trajectory predictionalgorithm of NEOs in order to be able to computeanalytically and numerically the minimum orbital intersection distances (MOIDs).The NEOPROP software consists of two separatemodules/tools: The Analytical Module makes use of analyticalalgorithms in order to rapidly assess the impact riskof a NEO. It is responsible for the preliminaryanalysis. Orbit Determination algorithms, as theGauss and the Linear Least Squares (LLS) methods,will determine the initial state (from MPCobservations), along with its uncertainty, and theMOID of the NEO (analytically). The Numerical Module makes use of numericalalgorithms in order to refine and to better assess theimpact probabilities. The initial state provided bythe orbit determination process will be used tonumerically propagate the trajectory. The numericalpropagation can be run in two modes: one faster(“ fast analysis ”), in order to get a fast evaluation ofthe trajectory and one more precise (“ completeanalysis ”) taking into consideration more detailedperturbation models. Moreover, a configurablenumber of Virtual Asteroids (VAs) will benumerically propagated in order to determine theEarth closest approach. This new “MOID”computation differs from the analytical one since ittakes into consideration the full dynamics of theproblem.
1. ARCHITECTURE
The general architecture of the tool, with special focuson the external and internal interfaces between thedifferent components, is assessed here.The high-level input/output interaction of the AnalyticalModule can be summarized as follows: As input, this module reads from an external file alist of observations; As output, it produces a file with the determinedinitial state and uncertainties of the NEO and thecomputed MOID;
As optional output, this module can produce a filecontaining a list of VAs, computed on the Line OfVariation (LOV).The high-level input/output interaction of the NumericalModule can be summarized as follows: As input, this module reads the main settings froman XML file (e,g, physical properties of the NEO,final epoch, integrator, etc). Moreover, it uses theoutput file created by the Analytical Module(containing the determined initial state) in order toinitialize the propagation; As optional input, it can import a list of VAs initialstates, precomputed by the Analytical module; As output, this module produces a file with thepropagated nominal trajectory and a separate outputfile with the list of close approaches to any planet. IfVAs trajectories are propagated, then an additionaloutput file provides their trajectory data; As optional output, it can produce a file with therelative positions of the NEO w.r.t. the third-bodiesused during the propagation. This file can beparticularly useful if the trajectory has to bevisualized.The described interfaces are summarized in Figure 1.NEOPROP Interfaces (note that the dashed lines areused for optional functionalities/files).
Figure 1. NEOPROP Interfaces
Having two stand-alone modules separates the two maintasks: orbit determination and trajectory propagation,allowing great flexibility. In this way, the differentfunctionalities can be combined in order to havedifferent analysis. Moreover, the integration of this toolwithin the NEO SSA Pilot Data Centre will bestraightforward (also thanks to the lack of dependenciesfrom external tools). This approach allows also to callthe Analytical and the Numerical Modules at the sametime for two different objects.Therefore, in order to perform the risk assessment of thefull NEOs catalogue (e.g. in the SBDC), only theAnalytical Module will have to be run for all theobjects. A list of PHAs will be determined, dependingon the value of the computed MOID (and absolutemagnitude). If the value of the MOID is critical (e.g. below 0.05 AU), the initial state will be numericallypropagated, along with VAs, in order to evaluate thedynamical effects on the MOID.When new observations will be available for a certainobject, only the analytical module will have to be runagain.In the next figure the high-level NEOPROP architectureis shown:
Figure 2. NEOPROP Architecture
A general aspect that should be highlighted is that boththe Analytical and the Numerical modules supportmulti-threading. The tool is able to detect how manythreads are available and to assign them in order toparallelize some computations. In the next paragraphsmore details are provided.
2. ANALYTICAL MODULE
The Analytical Module allows determining the initialstate (with uncertainty) of a NEO from a set ofobservations. First a Gauss algorithm is run in order to get apreliminary solution of the orbit determination problem.For this solution only optical observations are used. Ifthe test case has been run already once and a solutionalready exists, then the preliminary solution is notrecomputed, but it is imported (along with the alreadyused observations) and improved with the LLStechnique by means of the new set of observations.The Gauss method needs only three observations inorder to compute a preliminary initial state. The first,last and intermediate observations are by default theirst ones to be used. Then the tool runs the followingcombinations: The first and the last observations are kept fixed,while the third one is moved along the observationalarc; "forward loop": the last observation is used with anyother two; "backward loop“: the first observation is used withany other two.Before accepting any solution, the tool checks if thedetermined orbital elements are compatible with theorbit of a NEO (perihelion < 1.3 AU). The algorithmstops when a minimum RMS (predefined by the user) isreached or when all the combinations have beenexplored. In order to speed up this time-consumingprocess, a multi-threading approach has beenimplemented. In this way several combinations(depending on the available cores) are checked at thesame time. Figure 3. Gauss method algorithm
Then the preliminary solution is improved by means ofa LLS technique, which eventually considers also radarobservations. If requested, the initial states of a set ofVAs are computed along the LOV (two threads can beused for this computation). Finally, the determinedinitial state of the NEO is used in order to analyticallycompute its MOID by using the Sitarski method. The described work logic is showed by the followingdiagram:
Figure 4. Analytical Module detailed architecture
3. NUMERICAL MODULE
The Numerical module performs numericalintegrations. Its “fast analysis” is supposedto preliminary assess the trajectory of aNEO. When the MOID computed by theAnalytical module is small enough to bepotentially dangerous for the Earth,numerical integration is recommended.During this fast analysis, the initial stateprovided is numerically propagated and theclose approaches computed. Even though the “fast analysis” shares thesame models and XML input file with thoseused in the “complete analysis”, it does notallow the same flexibility. For instance, afixed set of perturbations is used during thisanalysis, while in the “complete analysis”the user can freely customize them. Thisapproach has been chosen for two reasons:1. The fast analysis is supposed to be“fast”. For this reason only some basic(and most important) perturbations areconsidered;2. The fast analysis may be run even byusers not so familiar with numericalpropagations. For this reason manysettings are initialized by default.The “complete analysis” is intended to beused for a complete risk assessment of aNEO. Starting from the orbit determined bythe Analytical module and from a set of VAs,the nominal trajectory will be propagatedand the “dynamical” MOID computed. VAstrajectories will be also propagated in orderto give to the user some more informationabout the sensitivity of the MOID due touncertainties in the initial conditions. Theinitial states of the VAs can be computed bythe Analytical Module along the LOV or by aMonte Carlo run. In this way a mix of randomVAs and VAs along the LOV can bepropagated.
In the Numerical module, the propagator reads the inputfrom XML file and initializes all the models(perturbations, integrator, asteroid, etc). Then thenumerical integration is started using as input the initialstate of an external file provided by the Analyticalmodule. During the propagation, all close approacheswith any planets are computed and stored. If theNumerical module is run in the "complete analysis"mode, then the VAs are initialized to be ready for thepropagation. Then, the propagator is called again foreach VA that needs to be propagated. If a multi-coreprocessor is used, VAs can be propagated in parallel,ince these computations are independent of each other.The work logic of the Numerical module is shown in thenext figure:
Figure 5. Numerical Module detailed architecture
The following external libraries/database are used inthis module: Spice ephemeris, asteroid orbits database(Astorb), time standards conversion tables (IERS), Earthorientation data (IERS), observatories databaseincluding their accuracy (MPC) and Earth sphericalharmonics (EGM96). In Table 1 and Table 2, the implemented integrators andperturbation models are listed.
Table 1. NEOPROP implemented integrators
Integrator Single/Multi-Step Step-SizeRunge-Kutta 45 single variable
Dormand Prince 8 single variable
Runge-Kutta 853 single variable
Runge-Kutta 4 single fixed
Runge-Kutta 4
Adapted single fixed*
Gauss-Jackson 8 multi fixed
Gauss-Jackson 8
Adapted multi fixed*
Gauss-Jackson 8
Self-Adapted multi fixed* *The integration follows a fixed step-size scheme, but for sometrajectory arcs (e.g. close to a celestial body) the step-size might bereduced by a factor of 10.
The term “adapted” has been used to distinguish theoriginal integrator scheme from a new algorithmimplemented in NEOPROP to improve the performance.The “adapted” integrators are fixed step-size integratorswhich are able to reduce the step-size whenever a certain condition applies. For the “adapted” integrators,the step-size is reduced by a factor of 10 when the NEOexits the sphere of influence of the Sun. This algorithmhas been implemented, because in these arcs theperturbing accelerations become significant and in orderto improve the accuracy a smaller step-size is required.On the other hand, along most of the trajectory (whengravity perturbations are not so relevant) larger step-sizegives better performance. The Gauss-Jackson 8 “self-adapted” integrator reduces the step-size when the ratiobetween the acceleration due to perturbations and thetotal acceleration is larger than 0.0001.Many physical properties can be specified tocharacterize the NEO: its mass, absolute magnitude,diameter, bulk and surface density, thermal conductivity,thermal capacity, albedo, infrared emissivity, spin axisobliquity and rotational period (no tumbling or binaryasteroid so far considered). At least the mass and thediameter or the absolute magnitude must be provided(the others can be set by default). These parametersaffect the solar radiation pressure perturbation, while theothers affect also the Yarkovsky computation. The list of implemented perturbations for the “fast” and“complete” analysis is summarized in Table 2.
Table 2. Implemented perturbation models “Fast” Analysis “Complete” Analysis
Third-Body (planets) Third-Body (planets + 4main asteroids)- Relativistic EffectsSolar Radiation Pressure Solar Radiation Pressure- YarkovskyEarth Spherical Harmonics Earth Spherical Harmonics
4. TEST CASES
The Analytical Module has been tested with 6 differentasteroids: the orbital elements and the MOID have beencomputed at certain reference epochs. As input theavailable observations (optical and radar) before thereference epoch has been used. As reference theNEODyS [2] system was used.
The 6 NEOs chosen for this test case are:2011AG5, 2012DA14, 2007VK184, 2009FJ,2004DC and (99942) Apophis. These havebeen selected because they all have a smallMOID (below 0.01 AU) and because theycover different types of orbits. Moreover, forthe NEOs 99942, 2004DC and 2012DA14also radar observations are available.For the asteroid (99942) Apophis 1681ptical observations and 7 radarmeasurements have been used.In order to have comparable results, inNEOPROP the output has been explicitlyrequested at the reference epoch 56400.0MJD (= 2456400.5 JD). This causes a smallloss in accuracy, since the initial state isdetermined at a certain epoch (middleobservation used in the Gauss method) andthen it is propagated up to the requestedepoch. The orbital fitting has been executed in twodifferent runs. This approach is also followedby the Orbfit tool, which is used in NEODyS.This type of approach is required when manyobservations are available and when theycover a very long observational arc. In thistest case we have 1688 observations spreadover a period of almost 10 years. The “multi-step fitting” strategy allows to initiallyconstrain the orbit using a smaller set ofobservations and then to further refine andimprove the determined initial state by usingthe full set of observations. For the NEOs 2011AG5, 2007VK184 and2009FJ, all the available observations (222,102 and 131) are used at once. Thepreliminary solution and orbital fit areperformed in a unique run. 2004DC requires a 2 steps orbital fit (164optical and 8 radar observations for the firstrun and other 35 optical for the second run).This was required due to the fact that these207 observations cover a very long arc(almost 8 years).2012DA14 requires a 3 steps orbital fit (191optical for the first run, 153 for the secondone and other 572 optical and 7 radarobservations for the last run). This approachis due to the close approach that this NEOhad with the Earth on the 15 th of February2013 [1]. This close approach changedconsistently its orbital elements. Moreover, alot of new observations were taken during it.For this reason the orbital fit processbecame quite complex. In this test case areference epoch before this close approachdate has been chosen. In this way the MOIDcomputation should be able to assess theclose approach distance reached on the 15 th of February.These last two NEOs together with Apophisallow simulating a re-computation of theorbital parameters due to new availableobservations. This is the typical scenario thetool will have to deal with. Only when a NEO is newly discovered, the total set ofobservations are used together (single ormultiple runs) to generate a first solution.Afterwards every time new observations aretaken, the solution is just updated by usingthe new measurements.In the next two tables, the results obtainedby NEOPROP are compared to the referencedata. Table 3. NEODyS-NEOPROP resultscomparison (99942, 2009FJ, 2011AG5)
Epoch [JD]
Δ-a [AU] %-a
Δ-e [-] %-e
Δ-i [°] -0.0016318 0.0012234 %-i
Δ-Ω [°] -0.0114062 0.0188884 %-Ω -0.006% -0.003% 0.014%
Δ-ω [°] %-ω
Δ-M [°] -0.0912437 -0.0000397 %-M
Δ-MOID [AU] %-MOID -0.001% 0.064% -1.040%
Table 4. NEODyS-NEOPROP resultscomparison (2012DA14, 2007VK184,2004DC)
Epoch [JD]
Δ-a [AU] -0.0000003 -0.0000100 0.0000000 %-a
Δ-e [-] -0.0000079 -0.0000031 0.0000000 %-e -0.007% -0.001% 0.000%
Δ-i [°] %-i
Δ-Ω [°] %-Ω
Δ-ω [°] -0.0042552 0.0785228 -0.0049840 %-ω -0.002% 0.107% -0.003%
Δ-M [°] %-M
Δ-MOID [AU] -0.0000342 0.0000012 -0.0000043 -MOID *computed with the Numerical module
For the asteroid 2012DA14, the valuereported as “MOID” is actually the closeapproach distance reached on the 15 th ofFebruary [5]. The MOID value found onNEODyS was not used, since it was based onthe new orbital elements determined afterthis Earth close approach (at epoch2456400.5) and not on the referenceelements used in this case (at epoch2456275.0, before the close approach date).This is why the closest Earth approachdistance was used as MOID.To check this value (which is not a MOID) theNumerical module was used, instead of theAnalytical one, to propagate the determinedinitial state until the day of the closeapproach. The computed numericalminimum close approach distance was equalto 0.00022899 AU, only 0.6% bigger thanthe reference value. The asteroid (99942) Apophis has been used to test theNumerical Module for a long-term propagation. Itsorbital elements (previously determined) have beenpropagated until 14-04-2029 in order to cover its nextEarth approach (on the 13-04-2029) and its consequentMoon close approach.
Only these two closeapproaches will be evaluated, sinceNEOPROP records only close approachesbelow 0.005 AU or below 100 times theradius of each celestial body. The twomentioned close approaches are the firstones which match this requirement and theyrequire about 25 years of numericalpropagation, which is enough to evaluatethe performance of the tool.
Since the initialstate used for the propagation comes from the Analyticalmodule, this test will allow evaluating the overallperformance of the NEOPROP software.The JPL Small-Body Database was used as reference forthe close approach distance [4], while the Horizons-JPLsystem was used as reference for the final state (asCartesian elements) [3].The final state and the close approaches distances (andepochs) are in line with the expected values. The resultscomparison is showed in the next table.
Table 5. JPL Horizons-NEOPROP resultscomparison (99942)
Δ-Value %-Value x [km] -2.139E+03 -4.366E+01 -0.02%Earth Distance [AU]
In the next table all the close approaches(from [4]) between 2006 and 2029 arereported, while in table Table 7 the dataobtained by NEOPROP is summarized andcompared to the reference values. Theresults of the close approaches are in linewith the expectations.
Table 6. Reference Apophis CloseApproaches [4]
Time (TDB) TimeUncertainty(hh:mm) Body Distance (AU)2006-Apr-10 23:49 < 00:01 Earth 0.202819844 < 00:01 Earth 0.096661113 < 00:01 Venus 0.078241899 < 00:01 Earth 0.216276061 < 00:01 Earth 0.11265166 < 00:01 Venus 0.124432355 < 00:01 Earth 0.000256194 < 00:01 Moon 0.000636262
Table 7. Apophis Close Approachescomputed by NEOPROP
Time (TDB) Body Distance(AU) Δ-Distance(AU) %-Distance(AU)2006-Apr-10,23:49
Earth 0.202820223 3.794E-07 0.000%
Earth 0.096665441 4.328E-06 0.004%
Venus 0.078240242 -1.656E-06 -0.002%
Earth 0.21627714 1.078E-06 0.000%
Earth 0.112651436 -2.238E-07 0.000%
Venus 0.124438646 6.291E-06 0.005%
Earth 0.000254828 1.366E-06 0.533%
Moon 0.000636175 8.743E-08 0.014%
This test case was also used to validate theimplemented integrators and to comparetheir performance. The last test case hasbeen rerun with each integrator using thesame set of inputs. Not all the integratorsettings are used by each integrator (e.g.the fixed step-size integrators ignore theminimum step-size info). These settings aresummarized here: able 8. Integrators settings
Parameter Value
Max step size 8640.0 sMin step size 0.01 sLocal Tolerance 1.0E-13Output interval 86400.0 s
The results have been compared with theoutput obtained in the previous test. Sincethe Gauss-Jackson 8 Adapted was used, thiswill be our reference. In the next table the results are collected.Only the Earth and Moon MOIDs have beencompared, since these represent the mostsensitive data. In this way it was alsopossible to reduce a bit the output to beshown.
Table 9. Integrators comparison
Δ-MOIDEarth %-MOIDEarth Δ-MOIDMoon %-MOIDMoon
RK45
DP8
RK853 -4.5716E-08 -0.02% 9.6784E-08 0.02%
RK4 -2.0973E-05 -8.26% -5.9186E-06 -0.93%
RK4_adapt
GJ8 -2.4058E-05 -9.47%
GJ8_S_adapt
All the integrators have very similar results,except for the fixed step-size integrators:RK4 and GJ8. Their worse performance isintrinsic to the nature of these integrators.Since they lack of any control on the localerror, they are not able to tackle situationswhere the perturbations become quite big(e.g. during a close approach with anothercelestial body). In order to improve theiraccuracy the fixed step-size should bereduced increasing the computational time.That was the reason why the “adapted”fixed step-size integrators have beendeveloped within this project and seem toperform better than their original integrationschemes.A last test case was run for the 2007VK184asteroid. The results of the orbitdetermination process have been alreadyshowed in paragraph 4.1. According to [6],2007VK184 has a possible Earth impactor in2048 (see details in the table below).
Table 10. 2007VK184 Earth Impactor data[6]
Time (TDB) σ (LOV) Distance (EarthRadius)2048-Jun-03 02:08
Therefore, the Analytical Module has beenrerun in order to generate VAs along theLOV. 2000 VAs have been generated and amaximum σ of 2 has been specified. Thenthe VAs has been passed to the Numericalmodule and propagated till 2048-Jun-04.Also NEOPROP detected an Earth impactorat exactly the same epoch. The minimumEarth distance is slightly different from thereference value; maybe due to differentvalues for the physical properties of the NEO(the ones used by NEODyS are unknown).Also the σ value of the VA is not the same;maybe due to a different definition of LOVadopted. Table 11. 2007VK184 Earth Impactor dataNEOPROP
Time (TDB) σ (LOV) Distance (EarthRadius)2048-Jun-03 02:08
5. CONCLUSIONS
The NEOPROP tool is able to determine the orbitalelements of a NEO starting from a set of observationsand to propagate them in order to compute the minimumclose approach distances to any planet. Particular focushas been put on determining impact risk with the Earth.The analytical MOID can be computed and used asmain driver for a preliminary analysis. If an advancedanalysis is requested, a set of VAs can be computed andpropagated in order to find out if an Earth impact mightoccur in the future.This software has been tested with several differentNEOs and all the results are within 1% of theirreference (NEODyS, JPL Small-Body Database andJPL-Horizons).
6. REFERENCES6. REFERENCES