Method and system of conjunction prediction between at least two bodies for collisions avoidance and rendez-vous manoeuvres

Abstract

The present invention is directed to a method (1000-6000; 1000A-E - 6000A-E) of predicting a conjunction between two bodies (50, 60; 50A,B; 50C,D; 50E,F, 50H). The method comprises the steps of: acquiring (2001; 2001A-G - 2001A-G) data defining at least one last known state and last known covariance of the two bodies; computing (2003-2007; 2003A-G - 2007A-G) a plurality of Gaussian Mixture Elements (GMEs) based on the acquired data; computing (3001-3003; 3001A-G - 3003A-G) a Gaussian Mixture Model (GMM) of an initial state based on the GMEs, the initial state being a state of the bodies at an initial time instant defining a lower boundary of a conjunction time interval comprising a time of closest approach of the two bodies; computing (3005-3015; 3005A-G - 3015A-G) Taylor expansions of dynamic equations of motions of the bodies for each GME in the GMM; computing (3017-3025; 3017A-G - 3025A-G) a final state and a GMM of the final state of the two bodies based on the Taylor expansions, the final state being a state of the bodies at a final time instant defining a higher boundary of the conjunction time interval; computing (3039; 3039A-G) a conjunction probability rate as a function of time based on the initial and final states and GMMs of the initial and final states, and computing (3045; 3045A-G) a conjunction probability by integrating the conjunction probability rate over the conjunction time interval, and determining (1015; 1015A-G) a possible conjunction if the conjunction probability exceeds a predetermined threshold. Advantageously, the step of computing (3017-3025; 3017A-G - 3025A-G) a final state and a GMM of the final state of the two bodies comprises: defining (3019; 3019A-G) a first matrix and a second matrix based on the Taylor expansion of state dynamic equations, the first matrix comprising coefficients of the Taylor expansion and the second matrix comprising exponents of the Taylor expansion, defining (3021; 3021A-G) a covariance matrix based on the GMM of the initial state, and computing (3023; 3023A-G) the GMM of the final state based on a combination of the first matrix, the second matrix and the covariance matrix.

Description

[0001] METHOD AND SYSTEM OF CONJUNCTION PREDICTION BETWEEN TWO OR MORE BODIES FOR COLLISION AVOIDANCE AND RENDEZ-VOUS MANOEUVRES

[0002] DESCRIPTION TECHNICAL FIELD

[0003] The present invention refers to the field of conjunction predictions of two or more bodies moving relatively from one to another. Particularly, the present invention relates to a versatile method and system of conjunction prediction that can be exploited in various fields of technology relating to moving bodies. For example, the present invention is referred to a method and system to predict and reach or avoid conjunctions between a spacecraft and another object, e.g., another spacecraft or space debris, respectively.

[0004] STATE OF THE ART

[0005] Conjunction prediction between two bodies is a useful information for any activity that involves moving an object, e.g. a ground vehicle, an aircraft, a spacecraft, etc. Indeed, the capability of predicting a possible conjunction in an effective and reliable manner allows preventing collision accidents both in manned and unmanned apparatuses. It should be noted that in the present description the term "conjunction” entails collisions, near misses, and rendezvous between two bodies in space, sky, land or water regardless of their type.

[0006] For example, more than nine hundred thousand artificial objects orbit the Earth and the probability of their conjunction with operational satellites keeps increasing with time, as new launches are performed, and in-orbit explosions and collisions keep happening. Managing close conjunctions in space and avoiding collisions is becoming a fundamental task for satellites and other spacecraft operators. Accordingly, there is a growing need in the space industry for automated and enhanced conjunction avoidance capabilities.

[0007] A generic conjunction avoidance process can be divided into three main steps:

[0008] 1. tracking of space objects;

[0009] 2. risk metrics calculation ( / .e., probability of conjunction and miss distance), and

[0010] 3. conjunction avoidance manoeuvre planning and execution.

[0011] These steps are partially executed by operators and / or systems that are based on ground, and partially by systems onboard the satellites.

[0012] Generally, the risk associated with a conjunction between two or more objects may be estimated according to the type of conjunction, each of which is generally identified as a short-term conjunction or as a long-term conjunction.

[0013] In detail, short-term, or high-velocity, conjunctions occur in Low Earth Orbit (LEO), where the relative velocity between the two objects is very high - e.g., up to 15 km / s. The known methods that address the evaluation of conjunction probability exploit a plurality of assumptions - allowed by the short-term features of above - to simplify calculating the conjunction probability.

[0014] Long-term or low-velocity conjunctions occur at higher orbits, e.g., Geostationary Earth Orbit (GEO), with respect to short-term conjunction and where the relative velocity is low {e.g., if both satellites are in very similar orbits). In this case, the assumptions made for short-term conjunctions do not apply and two main alternative approaches have been developed.

[0015] The first approach is based on a two-dimensional method, usually Patera's method due to the flexibility provided by using a line integral, as described in Patera, R. P. (2001). “General Method for Calculating Satellite Collision Probability1, published in Journal of Guidance, Control, and Dynamics Vol. 24, No. 4, pp. 716-722.

[0016] Other examples of method addressing short-term conjunctions comprise:

[0017] Foster, J. and Estes, H. (1992). "A Parametric Analysis of Orbital Debris Collision Probability and Manoeuvre Rate for Space Vehicles", NASA, National Aeronautics and Space Administration, Lyndon B. Johnson Space Center;

[0018] Patera, R. P. (2005). “Calculating Collision Probability for Arbitrary Space Vehicle Shapes via Numerical Quadrature", published in Journal of Guidance, Control, and Dynamics Vol. 28, No. 6, pp. 1326-1328;

[0019] Alfano, S. (2005). "A Numerical Implementation of Spherical Object Collision Probability", published in Journal of the Astronautical Sciences Vol. 53, No. 1, pp. 103-109, and

[0020] Chan, K. (1997). “Collision Probability Analyses for Earth Orbiting Satellites", published in 7th International Space Conference of Pacific-Basin Societies. Vol. 96, pp. 1033-1048.

[0021] For these cases, the conjunction geometry is defined by a set of volumes (cylinders, parallelepipeds, etc.) that can be divided into planar sections. In each time step, the two-dimensional probability is computed for each section, assuming a linear motion. Moreover, a one-dimensional probability is computed along the relative velocity vector on the volume. For each volume, the collision probability is obtained by multiplying the two probability values for the two volumes. Finally, the total collision probability is obtained by adding up the probabilities for the individual volumes.

[0022] A common limitation of these methods is that both objects must have a non-negligible relative velocity, which is a condition that is not always satisfied.

[0023] The second approach consists of computing a three-dimensional integral of the collision probability over time. This is done by defining the variable of collision probability rate and calculating the flow of this variable over the surface of hard-body volume during the conjunction time interval. Finally, the collision probability rate is integrated over time. Examples of methods addressing long-term conjunctions comprise:

[0024] Patera, R. P. (2003): “Satellite Collision Probability for Nonlinear Relative Motion" published in Journal of Guidance, Control, and Dynamics, Vol. 26, No. 5, pp. 728-733;

[0025] Alfano, S. (2007): “Beta Conjunction Analysis Toot’, published in AAS / AIAA Astrodynamics Specialist Conference. Vol.

[0026] 07, pp. 19-23;

[0027] McKinley, D. (2006): “Development of a Nonlinear Probability of Collision Tool for the Earth Observing System", published in AIAA / AAS Astrodynamics Specialist Conference and Exhibit. Eprint: https: / / arc.aiaa.Org / doi / pdf / 10.2514 / 6.2006-6295, and

[0028] Coppola, V. and McAdams, J. V. (2012): “Including Velocity Uncertainty in the Probability of Collision between Space Objects" published in AAS12247, 22nd, Spaceflight mechanics 2012, Vol. 143, pp. 2159-2178.

[0029] Due to various assumptions used to simplify the description of the orbiting bodies, the known methods to estimate collision probability are affected by a variety of limitations. All these methods reach a limited accuracy, particularly the collision probability is often estimated to a much higher value than the actual value thereof. The same applies to the estimation of the conjunction time between two objects. Both these values are key to correctly cope with an collision, possibly, by implementing a collision avoidance manoeuvre. In the art, a collision avoidance manoeuvre is recommended when the risk of collision is higher than a threshold commonly set 10’4. The small magnitude of the threshold leads to an immense number of collision warnings being generated by the known methods - e.g., up to hundreds of warnings per satellite per week -, while typically only a few collision avoidance manoeuvres per year are actually necessary. Thus, the known methods generate a very large number of false alarms and there is not any known solution that reliably solves this problem or is able to reliably assess whether a collision avoidance warning should require a follow-up action - e.g., collision avoidance manoeuvre - or not.

[0030] These drawbacks are further aggravated by the uncertainties in the orbit and propagation environment determination of bodies orbiting Earth. Uncertainty is introduced in the computations, from multiple sources. Primarily there is an uncertainty in the spacecraft state ( / .e., position and velocity), which is limited by the orbit determination process. Current techniques allow to determine the position of active satellites to centimetre accuracy. Conversely, space debris are more difficult to track and their position is determined with a precision in the order of tens to hundreds of meters. In the art, a plurality of uncertainty propagation methods has been proposed, which can be split into linear, nonlinear, and hybrid methods. Nonetheless, the Applicant has found that none of the solutions known in the art offers a fast and accurate method to foresee and prevent collisions between spacecrafts and other objects orbiting Earth.

[0031] SUMMARY OF THE INVENTION

[0032] It is an objective of the present invention to overcome the drawbacks of the prior art.

[0033] Particularly, it is an object of the present invention to provide a method and a related system adapted to provide an accurate and precise estimate of a conjunction probability between two or more objects and to implement collision or near misses avoidance manoeuvres, or rendez-vous manoeuvres accordingly. It is another object of the present invention to provide a method and a system for conjunction prediction and reaction adapted to reliably prevent undesired collisions or near misses and support planned rendez-vous both in case of short term conjunctions and long term conjunctions, in a fast and resources-effective manner.

[0034] It is a further object of the present invention to provide a method and a system adapted to provide an accurate and precise estimate of a reliable conjunction probability, particularly of spacecrafts having complex structures, in a fast and resource-effective manner.

[0035] These and further objects of the present invention will be clearer from the following description and from the annexed claims, which are an integral part of the present description.

[0036] According to a first aspect, the invention therefore relates to a method of predicting a conjunction between two bodies. The method comprises the steps of:

[0037] acquiring data defining at least one last known state and last known covariance of the two bodies;

[0038] computing a plurality of Gaussian Mixture Elements (GMEs) based on the acquired data;

[0039] computing a Gaussian Mixture Model (GMM) of an initial state based on the GMEs, the initial state being a state of the bodies at an initial time instant defining a lower boundary of a conjunction time interval comprising a time of closest approach of the two bodies;

[0040] computing Taylor expansions of dynamic equations of motions of the bodies for each GME in the GMM; computing a final state and a GMM of the final state of the two bodies based on the Taylor expansions, the final state being a state of the bodies at a final time instant defining a higher boundary of the conjunction time interval; computing a conjunction probability rate as a function of time based on the initial and final states and GMMs of the initial and final states, and

[0041] computing a conjunction probability by integrating the conjunction probability rate over the conjunction time interval, and

[0042] determining a possible conjunction if the conjunction probability exceeds a predetermined threshold.

[0043] Advantageously, the step of computing a final state and a GMM of the final state of the two bodies comprises: defining a first matrix and a second matrix based on the Taylor expansion of state dynamic equations, the first matrix comprising coefficients of the Taylor expansion and the second matrix comprising exponents of the Taylor expansion, defining a covariance matrix based on the GMM of the initial state, and

[0044] computing the GMM of the final state based on a combination of the first matrix, the second matrix and the covariance matrix.

[0045] The method according to the present invention yields faster, more precise and more accurate conjunction probability assessments with respect to the known methods, whilst keeping the computational load to a minimum. Further, the method according to the present invention can be used to estimate conjunction probabilities in any kind of conjunctions, particularly, in both long- and short-time conjunctions with no distinction. Particularly, by making use of matrices and matrix operations to compute the GMM of the final state, the computational load of determining the uncertainty propagation is largely reduced with respect to known methods.

[0046] In an embodiment of the invention, the step of defining a covariance matrix comprises:

[0047] defining a vector of initial state covariance values comprised in the GMM, and

[0048] defining the covariance matrix by compounding a plurality of vectors of the initial state covariance values, preferably the number of vectors comprised in the covariance matrix being equal to the maximum number of permutations of initial state covariance values.

[0049] In an embodiment of the invention, the step of computing the GMM of the final state based on a matrix combination comprises computing an expectation operator matrix as:

[0050] D — CO6ffSxC[ixn^comp]

[0051] where D is the expectation operator matrix, coeffs is the first matrix, ncompis the number of state variables, a C is a matrix defined as:

[0052]

[0053] where npermsis the maximum number of permutations of initial state covariance values, and B is a matrix defined as:

[0054]

[0055] I1 x

[0056] where A is the covariance matrix and exp is the second matrix.

[0057] In an embodiment of the invention, the step of computing the GMM of the final state based on a matrix combination comprises:

[0058] computing a covariance matrix P / j of the final state as:

[0059] ftft

[0060] pi ‘ - +p„<fc

[0061]

[0062] V?; Si" J

[0063] where Ci,pi...pnand Cj,qi...qnare coefficients of the Taylor expansion, p / y is the mean of the / -th and / -th components of a final state Xf, and 5xf1+(<1and 6xr”*cinare deviations, and E is the expectation operator, and

[0064] computing the mean , of the / -th component of a final state as:

[0065] ft =

[0066]

[0067] In an embodiment of the invention, the step of computing a plurality of GMEs comprises:

[0068] selecting a number of GMEs to be computed;

[0069] selecting a univariate splitting library based on the number of GMEs, and

[0070] scaling the univariate splitting library to fit a multivariate Gaussian distribution, defined by the initial states and by the initial state covariance of the two bodies.

[0071] In an embodiment of the invention, the step of computing a final state and a GMM of the final state of the two bodies comprises:

[0072] computing the Taylor expansion of the final state for a deviation from the initial state for each GME of the GMM, and computing a Taylor expansion of the final state for a deviation from the initial state for the GMM by combining the Taylor expansions computed for each GME.

[0073] In an embodiment of the invention, the step of computing a Taylor expansion of dynamic equations of motions of the bodies for each GME in the GMM comprises exploiting a Differential Algebra Computational Engine algorithm, which for each body is configured to:

[0074] defining the initial state of the body as a differential algebra variable in an inertial reference frame;

[0075] calculating environment variables as a function of the initial state;

[0076] defining aerodynamic forces, associated to the environment variables, in the inertial reference frame; computing an aerodynamic acceleration of the body based on the aerodynamic forces and properties of the body; computing gravitational acceleration due to Earth's gravity by determining respective Legendre polynomials and computing potential gradient in spherical coordinates;

[0077] computing a total acceleration by adding the previously computed accelerations, and

[0078] numerically integrating, by means of an integrator based on Runge-Kutta-Fehlberg (RKF).

[0079] In case of two bodies in space, the method comprises computing third body acceleration from Sun and Moon based on the positions of Sun and Moon retrieved from ephemeris in a SPICE library and use the third body acceleration together the other accelerations to compute the total acceleration.

[0080] In an embodiment of the invention, the step of computing a conjunction probability rate comprises selecting between a single sphere model and a multi-spheres model of each body, and performing a surface integration based on the selected model for each body.

[0081] Preferably, at least one multi-sphere model is computed by:

[0082] acquiring a three-dimensional model of the body,

[0083] fitting a plurality of spheres inside the three-dimensional body,

[0084] generating a sphere-mesh representation based on the plurality of spheres,

[0085] extracting a plurality of boundary points from the sphere-mesh representation, the boundary points defining a surface of the body, and

[0086] computing Lebedev weights, sphere centre positions and radii associated with the plurality of boundary points selected. A different aspect of the present invention regards a method of preventing a conjunction between two bodies in space. The method comprises:

[0087] executing the method of predicting a conjunction between two bodies in space according to any one of the previously described embodiments, and

[0088] if a possible conjunction is determined by the method of predicting a conjunction:

[0089] generating a corresponding conjunction data message;

[0090] computing a conjunction avoidance manoeuvre based on the conjunction data message, and command at least one of the two bodies in space to perform the conjunction avoidance manoeuvre.

[0091] The method according to the present invention substantially reduces a number of unnecessary conjunction avoidance manoeuvres performed thanks to the high precision of conjunction prediction attainable with the present invention. Moreover, conjunction data messages are provided in fast and timely manner, allowing defining an effective avoidance manoeuvre and, thus, ensuring the safety of a boy of interest, e.g.: a satellite, a space station, etc.

[0092] A further aspect of the present invention is directed to a system of conjunction prevention between bodies in space. The system comprises:

[0093] a conjunction prediction assembly, comprising at least a computing element, which is adapted to predict a conjunction probability between at least two bodies in space by implementing the method according to any of the previously described embodiments, and

[0094] a mission control assembly, comprising a computing element and communication element, which are adapted to compute and command to perform a conjunction avoidance manoeuvre to at least one body of the two bodies in space. In a different embodiment, the computing element and the communication element, of the mission control assembly, are adapted to compute and command to perform a rendez-vous manoeuvre to at least one body of the two bodies in space.

[0095] The system according to the above embodiment allows reaching a desired target via an optimized trajectory. For example, the method is used to determine the best trajectory for a spacecraft to reach a desired docking point, e.g. an orbital station or another spacecraft. Further, the method is used to perform controlled collisions, e.g. for disposal of decommissioned spacecrafts by reentry in the earth atmosphere.

[0096] In a different aspect of the present invention is directed to methods and systems for conjunction prediction between at least one aircraft, at least one ground vehicle, or at least one naval vessel and a second body (either artificial or natural) as defined in the attached claims. Particularly, the system comprises a conjunction prediction assembly, comprising at least a computing element, which is adapted to predict a conjunction probability between at least at least one aircraft, ground vehicle, or naval vessel, and the second body by implementing the method according to any of the previously-described embodiments, the system comprises a remote or an on-board mission control assembly, comprising a computing element and communication element, which are adapted to compute and communicate to the pilot of at least one of the aircrafts a collision avoidance manoeuvre to be performed, or command an autonomous or semi-autonomous control system of the aircraft for performing the collision avoidance manoeuvre.

[0097] Another aspect of the present invention is directed to methods and systems for conjunction prediction between at least one intercepting means - e.g., an artillery shell or a missile - and a static or moving target, as defined in the attached claims. Particularly, in the case of a static target the method comprises computing the conjunction probability as a Circular Error Probable (CEP) and the system is configured to command a gun laying adjustment of an artillery piece before firing the artillery shell or an adjustment manoeuvre of the missile in order to hit or avoid hitting the specified target.

[0098] Thanks to these embodiments, it is possible to greatly reduce the uncertainty in the area of impact of the intercepting means, or to prevent the intercepting means to hit an undesired target.

[0099] The systems according to the invention, grants the same advantages mentioned above with respect to the abovedescribed embodiments of the method mutatis mutandis.

[0100] Further characteristics and advantages of the present invention will be clearer from the following detailed description of some preferred embodiments thereof, with reference to the annexed drawings.

[0101] BRIEF DESCRIPTION OF THE DRAWINGS

[0102] The invention will be described here below with reference to some examples provided by way of example and not as a limitation and shown in the annexed drawings. These drawings show different aspects and embodiments of the present invention and, where appropriate, reference numerals showing like structures, components, materials and / or elements in different figures are denoted by like reference numerals.

[0103] Figure 1 is a schematic diagram of a system for conjunction prediction and collision avoidance according to an embodiment of the present invention;

[0104] Figure 2 is a flowchart of method of conjunction prediction and collision avoidance according to an embodiment of the present invention;

[0105] Figure 3 is a schematic block diagram of modules implemented by a conjunction prediction assembly according to an embodiment of the invention;

[0106] Figure 4 is a flowchart of a data precomputation method according to an embodiment of the invention;

[0107] Figure 5A to 5C are a flowchart of a conjunction prediction method according to an embodiment of the present invention;

[0108] Figure 6 is a flowchart of a Hafnian precomputing method according to an embodiment of the invention;

[0109] Figure 7 is a flowchart of a spacecraft modelling method according to an embodiment of the invention;

[0110] Figure 8 is a qualitative axonometric view of a three-dimensional spacecraft model used in the method of Figure 7; Figure 9 is a qualitative view of boundary points extracted from a multi-sphere model of a first portion of the spacecraft three-dimensional model of Figure 8;

[0111] Figure 10 is a qualitative view of boundary points extracted from a multi-sphere model of a second portion of the spacecraft three-dimensional model of Figure 8;

[0112] Figure 11 is a qualitative view of boundary points extracted from a multi-sphere model of the spacecraft obtained by combining the boundary points of the models of Figures 9 and 10;

[0113] Figures 12A and 12B show the conjunction probability rate and conjunction probability, respectively, output by the conjunction prediction method according to an embodiment of the invention for a first simulation, and

[0114] Figures 13A and 13B show the conjunction probability rate and conjunction probability, respectively, outputted by the conjunction prediction method according to an embodiment of the invention for a second simulation.

[0115] Figure 14 is a flowchart of method of conjunction prediction and rendez-vous execution according to an embodiment of the present invention;

[0116] Figures 15 and 16 are schematic diagrams of systems for conjunction prediction according to a second embodiment of the present invention, in which conjunctions between bodies in Earth's atmosphere - e.g., manned / unmanned aircrafts - are predicted;

[0117] Figures 17-22 are flowcharts of methods of conjunction prediction and avoidance / rendez-vous similar to those of Figures 2, 4, 5A-5C, 6, 7 and 14 for bodies in Earth's atmosphere according to the second embodiment of the present invention;

[0118] Figures 23 and 24 are schematic diagrams of systems for conjunction prediction according to a third embodiment of the present invention, in which conjunctions between bodies on land - e.g., manned / unmanned ground vehicles - are predicted;

[0119] Figures 25-30 are flowcharts of methods of conjunction prediction and collision avoidance / rendez-vous execution similar to those of Figures 2, 4, 5A-5C, 6, 7 and 14 for bodies on land according to the third embodiment of the present invention;

[0120] Figures 31 and 32 are schematic diagrams of systems for conjunction prediction according to a fourth embodiment of the present invention, in which conjunctions between bodies in water - e.g., manned / unmanned naval vessels - are predicted;

[0121] Figures 33-38 are flowcharts of methods of conjunction prediction and collision avoidance / rendez-vous execution similar to those of Figures 2, 4, 5A-5C, 6, 7 and 14 for bodies on water according to the fourth embodiment of the present invention;

[0122] Figure 39 is a schematic diagram of a system for conjunction prediction according to a fifth embodiment of the present invention, in which conjunctions between an artillery shell and a target are predicted;

[0123] Figure 40 is a schematic diagram of a system for conjunction prediction according to the fifth embodiment of the present invention, in which conjunctions a missile and a target are predicted;

[0124] Figures 41-48 are flowcharts of methods of conjunction prediction and target reaching similar to those of Figures 2, 4, 5A-5C, 6, 7 and 14 for an intercepting means according to the fifth embodiment of the present invention;

[0125] DETAILED DESCRIPTION OF THE INVENTION

[0126] Some preferred embodiments will be described in detail below, although the invention is susceptible to various alternative modifications. It must in any case be understood that there is no intention to limit the invention to the specific embodiment illustrated, but, on the contrary, the invention intends to cover all possible alternative or equivalent methods' steps or systems' devices that fall within the scope of the invention as defined by the attached claims. Unless otherwise defined, all the terms of the art, notations and other scientific terms used herein are intended to have the meanings commonly understood by those skilled in the art to which this description belongs. In some cases, terms with commonly understood meanings are defined herein for clarity's sake and / or ready reference; the insertion of such definitions in the present description must therefore not be interpreted as representative of a substantial difference with respect to what is generally understood in the art.

[0127] Particularly, the term 'conjunction' in the present description comprises collisions, near misses, and rendezvous between two bodies regardless of their types - e.g., spacecrafts, aircrafts, ground vehicles, naval vessels, space debris, static objects (either natural or artificial).

[0128] The terms "comprising”, "having”, "including” and "containing” are to be understood as open terms ( / .e. the meaning "comprising, but not limited to”) and are to be considered as a support also for terms such as "essentially consist of', "essentially consisting of”, "to consist of' or "consisting of”.

[0129] The use of "for example”, “etc", "or” denotes non-exclusive alternatives without limitation, unless otherwise noted. The use of "includes” or "comprises” means "includes / comprises, but not is limited to”, unless otherwise noted.

[0130] With reference to the block diagram of Figure 1, a system for conjunction prediction and collision avoidance according to an embodiment of the present invention, from now on simply identified as "system 1”, comprises a conjunction prediction assembly 10, connected to one or more data repositories, only one data repository 20 being shown in the example, and to a mission control assembly 30. Further a monitoring assembly 40 is connected to the data repository 20.

[0131] In detail, the conjunction prediction assembly 10 generally comprises one or more processors {e.g., CPUs, GPUs, ASICs, DSPs, etc.), volatile / non-volatile memory modules, communication modules {e.g, modems, routers, etc.), input / output interfaces {e.g., user interfaces, communication ports, etc.), and ancillary modules {e.g., power modules). As will be clear to the skilled person, the conjunction prediction assembly 10 can be implemented by a single machine or by a computer network comprising two or more physical and / or virtual machines.

[0132] The data repository 20 generally comprises one or more, physical and / or virtual, data storing unit that are configured to store large amount of information - typically as digital data, e.g. binary data, and to provide such information to one or more user, such as the conjunction prediction assembly 10. In the considered example, the data repository 20 stores space debris data, spacecraft data and operational data generated at least by the conjunction prediction assembly 10. The mission control assembly 30, as known, comprises ground-space communication apparatuses and computer networks that allows exchanging with a controlled spacecraft 50. Particularly, the mission control assembly 30 receives telemetry data regarding the spacecraft 50 and is adapted to provide operating instructions to the spacecraft 50 - e.g., operating instructions necessary to perform a desired manoeuvre. The telemetry data are, preferably, stored in the data repository 20 among the spacecraft data by the conjunction prediction assembly 10 or directly by the mission control assembly 30.

[0133] The monitoring assembly 40 of Figure 1 comprises ground-based radars and optical telescopes, space-based sensors, computer networks to track objects orbiting Earth and their orbits, such as the monitoring of spacecraft / space debris performed by European Space Surveillance and Tracking (EU SST) network, U.S. Space Command (USspaceCOM) and NASA Orbital Debris Program Office (ODPO). Preferably, the data generated by the monitoring assembly 40 are stored in the data repository 20 as the space debris data already mentioned. The space debris data may be updated asynchronously - e.g., as soon as one or more data values changes - or periodically with a predetermined time period. The spacecraft 50 comprises any vehicle / equipment deployed in orbit around the Earth or travelling through a predetermined region of space - e.g., from Earth's Thermosphere outwards. In a non-limiting embodiment of the invention, for each monitored spacecraft 50 by the system 1, the data repository 20 stores a three-dimensional (3D) model thereof, which will be described in detail later in the description.

[0134] Space debris 60, comprise any (artificial) object in orbit about the Earth that no longer serves any useful purpose -e.g., derelict spacecraft and upper stages of launch vehicles or parts thereof, debris intentionally released, debris created in explosions or collisions, solid rocket motor effluents, and material (such as flecks of paint) released from a spacecraft due to thermal stress or small particle impacts.

[0135] The system 1 executes a conjunction prediction and collision avoidance method 1000 of which Figure 2 is a flowchart, which estimates a probability between a spacecraft 50 and another object, such as space debris 60 and enforces at least one conjunction-avoidance manoeuvre to avoid conjunctions, i.e., at least one collision-avoidance manoeuvre is enforced to avoid a collision between the spacecraft 50 and the space debris 60.

[0136] The method 1000 is initiated by the conjunction prediction assembly 10 when a data update is received and / or a start command is received (e.g., generated by a user - either a human or an automated agent - or generated on a predetermined time-base) (input block 1001).

[0137] The conjunction prediction assembly 10 acquires or updates the debris data, the spacecraft data stored in the data repository 20 (step 1003). Alternatively, debris data and / or spacecraft data are directly acquired from the monitoring assembly 40 and from the mission control assembly 30.

[0138] Preferably, one or more data filters (and / or sieves) are applied to the debris data, in order to exclude data related to space debris that are not expected to intersect the orbit of the spacecraft 50. For example, one filter is based on a relative distance between the spacecraft 50 and the space debris 60 - e.g., if the spacecraft 50 moves in a substantially circular orbit at a first altitude (such as 800km) data about space debris 60 moving in a circular orbit at an altitude of 1000 km will be ignored since a conjunction between the two is highly improbable. In conclusion, the settings of data filters and / or sieves depend on the orbits radii and eccentricities or travelling path of the spacecraft 50 and of the space debris 60. At the end of this filtering phase, one or more pair spacecraft 50 - space debris 60 are determined (step 1005) and their conjunction probability is computed as disclosed in the following.

[0139] The conjunction prediction assembly 10 computes a conjunction probability for each space debris 60 and spacecraft 50 pair (step 1007).

[0140] The conjunction prediction assembly 10 determines if one or more of the computed conjunction probabilities exceed a conjunction threshold (decision step 1009). For example, the conjunction threshold is equal to a conjunction probability of 10-4.

[0141] If there is not any conjunction probability exceeding the conjunction threshold (output branch N of step 1009), the procedure is interrupted, since there is not any risk of conjunction. Accordingly, the conjunction prediction assembly 10 enters an idle state, waiting for a data update and / or a start command to reiterate the method from step 1001. If one or more of the conjunction probabilities exceed the conjunction threshold (output branch Y of step 1009), for each conjunction probability exceeding the conjunction threshold a corresponding conjunction data message is generated (step 1011). For example, a generic conjunction data message comprises information identifying and / or defining the spacecraft 50 and the space debris 60 associated with the warning, telemetry information - e.g., Time of closest approach (TCA), conjunction data message (CDM) creation date, trajectory, relative speed, etc. of the spacecraft 50 and of the space debris 60.

[0142] The conjunction data messages are provided to the mission control assembly 30 (step 1013), which in its turn analyses each conjunction data message to determine if a conjunction-avoidance manouevre is necessary (decision step 1015). In the negative case (output branch N of step 1015), any manoeuvre is not (yet) needed, and the method comprises waiting for a data update and / or a start command to reiterate the method from step 1001.

[0143] On the contrary, if a collision-avoidance manoeuvre is required (output branch Y of step 1015), the mission control assembly 30 determines the parameters characterizing the collision-avoidance manoeuvre (step 1017) to be performed by the spacecraft 50 to avoid the space debris 60. Then, the mission control assembly 30 send manoeuvre instructions to the spacecraft 50 to perform the collision-avoidance manoeuvre (step 1019) and monitors the execution thereof to its completion (step 1021). For example, the mission control assembly 30 stores one or more instructions sets that allows implementing one or more procedures to assess the necessity of collision-avoidance manoeuvre, to compute the collision-avoidance manoeuvre based on the conjunction data messages and to command the spacecraft 50 to execute the collision-avoidance manoeuvre.

[0144] Afterwards, the method 1000 proceeds to wait for a data update and / or a start command to reiterate the method from step 1001.

[0145] In a first embodiment of the invention, the conjunction prediction assembly 10 implements a plurality of computing modules, each of which is configured to perform computation that allows computing a conjunction probability between two bodies A and B - e.g., the spacecraft 50 and the space debris 60 -, in a precise, accurate and fast manner. In the exemplary blocks diagram of Figure 3 the conjunction prediction assembly 10 comprises the following modules: a Gaussian Mixture Model (GMM) module, or GMM module 101 in brief, a Differential Algebra (DA) module 103, an uncertainty propagation module 105 and a conjunction probability module 107. In an alternative embodiment of the invention, a single module, i.e. a DA-propagation module is provided, which replaces the DA module 103 and the propagation module 105 by combining the operations performed by the modules 103 and 105 as described below. In the considered embodiment of the invention, the conjunction probability between two bodies A and B is computed based on a hybrid Differential Algebra - Gaussian Mixture Model (DA-GMM) of the interaction between the bodies A and B. In this model, due to intrinsic uncertainty in the measurement of position and uncertainty affecting the trajectories computation each body A and B is modelled by a probability density function (PDF). In particular, the probability distribution describing the bodies A and B is defined by a Gaussian Mixture Model (GMM) that is a weighted sum of N Gaussian distributions, named Gaussian Mixture Elements (GMEs). Further, a final state - e.g., at a time of closest approach (TCA) between bodies A and B - of each body A and B is computed by using a Taylor series expansion of the dynamical equations of motion of the body - which are referred to also as ‘dynamics' for short. Advantageously, the final state is evaluated by computing a deviation from an initial state with Taylor expansion coefficients, instead of using floating point variables. Afterwards, the uncertainty propagation of each individual GME to the final state is performed. Finally, the conjunction probability is derived from the propagated probability density function.

[0146] Particularly, as used in the description, the term "state” indicates a set of information on the position in space of a body and its velocity - e.g., a state comprises position and velocity data, for example both expressed in set of three components referred to a three-dimensional reference system or reference frame:

[0147]

[0148] ,v ~ [ j- r !?; x y s sy ty f1

[0149] In details, the computing modules 101 - 107 execute a method 2000 of data precomputation and a method 3000 of conjunction prediction, described with reference to the flowchart of Figure 4, and 5A and 5B, respectively, in order to determine the conjunction probability between the bodies A and B.

[0150] The data precomputation method 2000, for example, is performed at least once, particularly when the conjunction prediction assembly 10 starts monitoring a new spacecraft 50. In general, input data and other variables used are stored by the conjunction prediction assembly 10, in the data repository 20 (and accessed by the conjunction prediction assembly 10) and / or input by a human or automatic software agent I / O interfaces.

[0151] In detail, the GMM module 101 acquires input data (step 2001) comprising:

[0152] - a number of Gaussian Mixture Elements (GMEs) to be used,

[0153] - most recent, or last, known state of the bodies A and B XA(IO) and XB(IO), where to is a most recent time instant to which tracking data used for the computation refers, and

[0154] - a most recent, or last, known state covariance of the bodies A and B PA(IO) and Rs(to), also referred to the most recent time instant.

[0155] It should be noted that each body A and B may have a respective most recent known time instant since the bodies may be tracked by different systems and or with a different time period.

[0156] Based on the number of GMEs a univariate splitting library is selected (step 2003), among a set stored in the memory of the conjunction prediction assembly 10 or the repository 20. The univariate splitting library comprises a Gaussian sum approximation of the unitary one-dimensional Gaussian distribution:

[0157]

[0158] where N is the number of Gaussian Kernels - i.e., the GMEs - anda, waand a2are the mean, weight and variance of the cr-th Gaussian density function pg.

[0159] For example, the univariate splitting is defined based on orwood, J., Aragon, N., and Poore, A. (2011): “Gaussian Sum Filters for Space Surveillance: Theory and Simulations" published in Journal of Guidance, Control, and Dynamics Vol. 34, No. 6, pp. 1839-1851.

[0160] Then, the univariate splitting library is scaled to fit a multivariate Gaussian distribution, pg(x,[j,P, defined by the initial states XA(IO), XB(IO) and by the initial state covariance PA(to), Pe(to) of the bodies A and B, respectively, (step 2005). Finally, the multivariate Gaussian distribution is stored in the data repository 20 and / or directly by the conjunction prediction assembly 10 (step 2007).

[0161] Preferably, the data precomputation method 2000 is performed at least once for each spacecraft 50 monitored by the system 1 and / or any relevant space debris 60 detected.

[0162] The conjunction prediction method 3000 (Figures 5A-5C) is then performed as follows. The GMM module 101 acquires updated input data (3001) comprising:

[0163] - the multivariate Gaussian distribution computed for bodies A and B by means of the method 2000, - an updated last known state of the bodies A and B XA(IO) and XB(IO), where to is last known time instant of propagation, and

[0164] - an updated last known state covariance of the bodies A and B PA(IO) and Ps(to).

[0165] Based on the updated input data, the multivariate Gaussian distribution is refined over different dimensions (steps 3003). In an embodiment of the present invention, the refinement over different dimensions refers to an extension from a univariate distribution to the multivariate distribution under consideration, as explained in Norwood, J. et al. mentioned above.

[0166] In summary, the GMM module 101 provides as output a GMM of an initial state with respect to the bodies A and B - i-e., {Pi, jJi, w,)to for the generic / -th GME of the GMM. The initial state is a state of the bodies at an initial time instant, or initial time to in the following, defining a lower boundary of a conjunction time interval comprising a time of closest approach (TCA) of the two bodies.

[0167] Preferably, the GMM module 101 acquires as input parameters an estimated time length of the conjunction - varying based on the orbital scenario, e.g. high or low velocity conjunctions - and an estimated TCA provided, for example by the monitoring assembly 40. The lower boundary of the conjunction time interval is then the first epoch of a conjunction time window surrounding the given TCA, whose length is, for example, predetermined. The last known states of the bodies, as well as the associated last known covariances, are propagated until the start of the conjunction, i.e. lower boundary of a conjunction time interval. Additionally, even though not limitatively, an updated TCA is computed by the GMM module 101 to be used instead of the TCA received as an input.

[0168] In the embodiment considered, the following steps are iterated for each generic time instant f comprised between the initial time to and a final time instant, or final time tffor short. The final time tf is the upper boundary of the conjunction time interval. Particularly, each time instant f is separated by a time step At from previous and / or following time instants t-i and tj+1. The time step At may be defined as having a substantially constant value in simpler embodiments of the invention. On the contrary, in embodiments having a higher complexity, the length of the time step At varies dynamically in time as a function of changes in the velocities of the bodies during the conjunction.

[0169] Accordingly, a time instant variable / it is initialized at the initial time to.enc of the conjunction time interval (step 3005). The Differential Algebra, DA, module 103 is adapted to compute a trajectory estimation of the bodies A and B up to a the time instant To this extent, the DA module 103 performs an integration of Taylor series in differential algebra. Particularly, the DA module 103 acquires (step 3007):

[0170] - a dynamics model coded in a DA software with analytical expressions,

[0171] - an order k of the T aylor expansion,

[0172] - the initial time to, enc, a final time tf and the time step At of a conjunction time interval - also referred to as integration time or integration period -, and

[0173] - properties of the bodies A and B - e.g., frontal areas AA, AB, masses mA, me, drag forces acting on the bodies CD, A, CD.B, etc.

[0174] Preferably, the conjunction time interval is defined by specifying a threshold value. For example, in 'regular' scenarios, where both bodies A and B have different orbits, a substantially punctual conjunction in time is detected. In the considered example, the conjunction probability rate monotonously increases until reaching the maximum value(s) at defined time instant(s). These maxima are associated with the areas in which the bodies A and B are at relative shortest distances - i.e., areas in which the bodies A and B are at the closest distances one another reached on their current trajectories. Away from these areas the distance between bodies A and B grows monotonously, while the conjunction probability decreases monotonously. Accordingly, the conjunction time interval is selected by setting a threshold value for a minimum conjunction probability rate that should be considered. Preferably, a threshold value equal to, or lower than, IO10, or more preferably a threshold value equal to 10-12. Alternatively, the conjunction time interval is selected based on a threshold value to ensure a nominal separation distance between the bodies A and B. This allows pre-computing the conjunction time interval before performing the simulations; therefore, reducing the computational load required to determine the conjunction probability. Instead, in different scenarios from the regular scenario - e.g., where the conjunction risk is not continuous -, the conjunction time interval cannot be simply selected by means of a probability / separation threshold as in the previous examples. On the contrary, the selection of the conjunction time interval is based on multiple (case-dependent) parameters, such as tracking station passes and available computational load. Particularly, a tracking station pass is any interval of time during which the relevant body - e.g., the spacecraft 50 - is in a direct line of sight of a tracking station - e.g., comprised in the mission control assembly 30. As should be clear to the skilled person, the conjunction time interval and / or the threshold value may be pre-selected or pre-computed for each couple of bodies considered, before implementing the method 3000, e.g. in a similar way as described in method 2000 for the multivariate Gaussian distribution.

[0175] Afterwards, the DA module 103 initialises a plurality of variables of a DA environment (step 3009) and an environmental and force model for the initial time of propagation (step 3011) based on the Taylor order k, the conjunction time interval and bodies A and B properties. The trajectory of the bodies is numerically integrated over a single time step At of the conjunction time interval in the DA environment (step 3013). The DA module 103, thus, computes, and preferably outputs, a Taylor expansion of the final state for a deviation from the initial state - due to uncertainties - for each / -th GME for the bodies A and B, e.g.

[0176]

[0177] and J’(XB* (t)) and, preferably, states of the bodies A and B as a function of time XA< (t) and xei (t), where t is an / -th time instant comprised between the initial time instant to.enc and a final time instant tf of the conjunction time interval. Finally, the DA module 103 computes the state XA(t) and Xs(ti) and a Taylor series expansion of the state of the bodies A and B at time t, with respect to the initial time to.enc for the whole GMM (step 3015). For example, the data regarding the initial states of the bodies A and B XA(IO), XB(IO) are initialised as DA variables: [XA,O] = XA(IO) +6xA(to), [XB,O] = XB(IO) + 5xB(to), where 5XA,B is a deviation - i.e. , a random vector - accounting for the uncertainty associated with the states of the bodies A and B.

[0178] The environmental and force model is coded in a differential algebra software, using analytical expressions, and is initialized based on the properties of the bodies A and B. Advantageously, the environmental and force model provides models describing the motion of the bodies, perturbing accelerations affecting the bodies, the gravity field acting on the bodies, shape and rotation of the bodies, third body effect and the atmosphere in which the bodies are immersed. Similarly, the Taylor expansions of the propagation to the final state are computed by a differential algebra method, preferably implemented by one or more algorithms based on the Differential Algebra Computational Engine (DACE) -i.e., a differential algebra library developed in C++ and created by Dinamica srl. for the European Space Agency (ESA), which is described in Rasotto, M., Morselli, A., Wittig, A., Massari, M., Di Lizia, P., Armellin, R., Valles, C., and Ortega, G. (2016): “Differential algebra space toolbox for nonlinear uncertainty propagation in space dynamics", published in 6thInternational Conference on Astrodynamics Tools and Techniques (ICATT).

[0179] In an embodiment, the DA module 103 is configured to compute an aerodynamic acceleration of a body, e.g. body A or the spacecraft 50. Firstly, the initial state of the satellite is provided as input and used as a differential algebra variable. This state may be received in one of different reference frames (e.g., J2000, ECEF or TEME) and is converted in an inertial reference frame such as J2000. The state is used to calculate secondary parameters (e.g., orientation, angles), which are in their turn, exploited to compute environmental variables (e.g., longitude and latitude, the flightpath angle and the heading angle of the body). These variables are used to determine the atmospheric density at a given location in space, i.e. in the bodies positions, and to transform the associated aerodynamic forces from the aerodynamic reference frame to the inertial reference frame. The transformed aerodynamic forces together with body properties (e.g., mass) are used to compute the aerodynamic acceleration. The gravitational acceleration due to the Earth - considering the spherical harmonics up to degree and order six - is computed analytically by determining the Legendre polynomials and computing potential gradient in spherical coordinates. The spherical gradient is converted to Cartesian coordinates to obtain the spherical harmonics acceleration. The third body acceleration from the Sun and the Moon is calculated once the positions of these bodies are retrieved from ephemeris in the SPICE library (Spacecraft Planet Instrument C-matrix Events). The accelerations are added together in a total acceleration, which is numerically integrated (e.g., by means of an integrator based on Runge-Kutta-Fehlberg (RKF) 7(8)) on a time step At to obtain the state XAj(ti) and x^(f) and the Taylor series expansion (xA)(tl)) and (XB)(t) of each state of the body A or B at time t with respect to the preceding time instant

[0180] Then, operation passes to the uncertainty propagation module 105, which is adapted to determine a Gaussian Mixture Model of the final state uncertainty (R, pi,Wi)tf.

[0181] Particularly, the propagated covariance matrix Pjj and the mean pi of the probability distribution of a random variable x, described by using the Taylor expansion, can be defined by an expectation operator E[.] of the random variable x and of its deviation 5x, respectively. Preferably, an eigenvalue clipping procedure is applied to the covariance matrix Pi, to ensure a positive semi-definite covariance matrix. Based on Valli, M., Armellin, R., Di Lizia, P., and Lavagna, M. R. (2013): “Nonlinear Mapping of Uncertainties in Celestial Mechanics" , published in Journal of Guidance, Control, and Dynamics Vol. 36, No. 1, pp. 48-63, the propagated covariance matrix P / j of a propagated state is calculated as:

[0182] PiPj

[0183]

[0184] where Ci,pi...pnand Cj,qi...qnare the Taylor coefficients of the Taylor polynomial expansion, fjij is the mean of the / -th and / -th components of a final state Xf, and 5xf1+(<1and bxp^ are deviations.

[0185] Further, the mean , of the / -th component of a final state y, is defined as:

[0186]

[0187] Solving equations (2) and (3), thus, involves computing the expectation E[.] of a product of normal random variables ( / .e., 5xiP1+c>1and 5xpn+cin, and 5xf1and 5xpn). Based on the teachings of Isserlis, L. (1918): “On a Formula for the Product-Moment Coefficient of any Order of a Normal Frequency Distribution in any Number of Variables" published in Biometrika Vol. 12, No. 1 / 2, pp. 134-139, the expectation is defined as:

[0188] . , I 0, if s is odd

[0189] h Ix'dx? ...x ,, f - <

[0190]

[0191] 1‘ | HaffPj, if sis even (4)

[0192] where s = si + ... + snis the order of each T aylor series expansion terms, and Haf(.) is the Hafnian operator. Considering that the equation (4) stems from a Taylor series expansion, the maximum value of s, i.e. sn, coincides with the Taylor expansion order k (sn= k). When s is even, the expectation of that term is equal to the value of the Hafnian operator evaluated for the initial state covariance matrix P. For this calculation, an expanded state vector is defined from the initial deviation 6x, as:

[0193] z- [-Si ...... ifils_, |

[0194]

[0195] where l's / is a unit vector having a size equal to the exponent of the corresponding deviation 5 , i.e. s,, and zj, Z2, Z3 are each a component of the deviation 5x, (i.e., the deviation from the initial state). The initial state covariance matrix P is associated with the expanded state vector and the Hafnian operator by a given matrix Z defined by components Oij is calculated as (based on the above-mentioned paper by Valli et al., 2013):

[0196] = £ [J

[0197]

[0198] . , (6)

[0199] where D, is the set of permutations of 1,2,..., s that satisfy the property pi < p3 < ps < ■ ■ ■ < ps-i and pi < p2, ps < P4,..., ps-i < ps. The complexity of the Hafnian operator is bound to the high number of function evaluations that are required as the number of permutations increases. As known in the art, the DA-GMM approach comprise computing the propagation of covariance for every Gaussian mixture element at each time step. Thus, DA-GMM procedures known in the art require a large computational load that becomes quickly unsustainable in a practical application (e.g., up to 100,000 evaluations per time step for a given conjunction scenario).

[0200] For example, let us assume a simplified scenario with bodies having a two-dimensional motion. The covariance matrix is 4 x 4 and for any given term of the summation in equation (2), the evaluation of the expectation operator E[.] results in the expression f(5x,P)

[0201]

[0202] _(7)

[0203] This operation can be performed for every summation term in the covariance propagation equation, resulting in a total of n2COeff expressions, where ncoeff is the number of terms in the T aylor series expansion of each component. The number of addends in these expressions will be equal to, or lower than, the maximum number of permutations of covariance values Pij, i.e. nperms.

[0204] In order to determine the GMM final state uncertainty in a fast and effective manner, the method 3000 comprises having the uncertainty propagation module 105 acquiring (step 3017):

[0205] - the Taylor expansions T(xA(ti)) and T(xB(ti)) computed by the DA module 103, and

[0206] - the initial state covariance of the bodies A and B PA(IO) and Ps(to).

[0207] Then, the uncertainty propagation module 105 defines precomputed Hafnian operator matrices (step 3019) based on the order k of the Taylor expansions T(xA(ti)) and T(xB(ti)) and the populate the Hafnian operator matrices based on the Taylor expansions T(xA(ti)) and T(xB(ti)) (step 3021). Particularly, the Applicant noted the values of the initial state covariance matrix changes for each uncertainty propagation to be computed, the Taylor expansion structure is constant. Based on this observation, the Applicant developed a Hafnian precomputing method 4000 which implements above-mentioned steps 2019 and 2021 - of which Figure 6 is a flowchart.

[0208] Taylor expansions data are acquired (step 4001). Preferably, the Taylor expansions (xA(ti)) and (xB(ti)) and their order / rare acquired.

[0209] Based on the order k of the Taylor expansions (xA(ti)) and (xB(ti)) exponents and coefficients of the Taylor expansions (xA(ti)) and (xB(ti)) are identified (step 4003).

[0210] Then, a exponents matrix and a coefficients matrix are created based on the identified exponents and coefficients (steps 4005 and 4007). The exponents matrix has n2coeff rows and npefms-n^omp columns, where ncompis the number of state variables XA(t) and XB^), six in the considered example. The coefficients matrix has n2coeff rows and npefmscolumns, where the value of each column, n represents the coefficient that multiplies the n-th addition term in equation (7). The exponents matrix comprises the exponent of the Taylor expansions and the coefficients matrix comprises the coefficients of the Taylor expansions (xA(ti)) and (xB(ti)) of each GME of the GMM.

[0211] By considering the two-dimensional example of equation (7) of above, the first rows of the coefficient and exponent matrices are:

[0212] coeffs =

[0345] (8)

[0213] exp = [1 1 00002 1 3 1 0 1], (9)

[0214] In summary, the method 4000 comprises defining two matrices, a coefficient matrix coeffs, with a size n2COeffxnPerms, comprising the coefficients of the Hafnian operator Haf(P) and an exponents matrix exp, with a size of n2COeffxriperms n2COmp, comprising the exponents of the Hafnian operator Haf(P). For each matrix, each column corresponds to an addition term in the summation operator of the covariance matrix Py (Equation (2)).

[0215] Advantageously, the coefficient matrix coeffs and the exponents matrix exp must be defined - or loaded from a memory or the repository 20 - only once by the uncertainty propagation module 105, thus, saving computational load and time while computing results for all the GME of the GMM over each time step comprised in the conjunction time interval. Further, thanks to this method, the full expectation operator and the summation along all Taylor expansion coefficients are computed as a matrices products - i.e., the complexity of computing the GMM final state uncertainty is simplified being substantially limited to matrix operations.

[0216] By making use of the precomputed matrices of the Hafnian operator, as described above, for a given expansion order and number of components in the coefficients and exponents matrices, the computational load of determining the uncertainty propagation is largely reduced, as apparent from the test data in Table 1 herein below. Tests performed by the Applicant show that the computation of the uncertainty propagation is more than 3,700 times faster using the method herein described w.r.t. a full Hafnian computation known in the art, for Taylor expansion order k = 3. Furthermore, the computational speed is more than 22,000 times faster using the method herein described w.r.t. a full Hafnian computation known in the art, for Taylor expansion order k = 4.

[0217] Taylor expansion order Full Hafnian Matrix simplification

[0218] 190 s 0.051 s

[0219]

[0220] 4 >10 h 1.63 s

[0221] Table 1

[0222] Thanks to this massive reduction in computational time, it is possible to estimate the conjunction probability in a timeeffective manner from simple models - which include a limited number of GMEs and a short conjunction time - to complex models - which include a large number of GMES and long conjunction periods.

[0223] Back to the method 3000, subsequently to computing the Hafnian operator for the Taylor expansions (xA(ti)) and (xB(ti)), the uncertainty propagation module 105 computes the final state uncertainty for each / -th GME of the GMM (step 3023) and, accordingly, computes the final state uncertainty for the GMM (step 3025). The process followed to reconstruct the mean and covariance of the propagated final state is based on equations (2) and (3). Particularly, the initial state covariance P,j is redefined as a vector Avecrepeated npeimstimes, which in the considered two-dimensional example is expressed as:

[0224] \PUP12 P21 P22I

[0225]

[0226] The expectation operator E[.] (equation (4)) is computed with matrix operations:

[0227] B = A™? 2n'iper ms ,l

[0228] f'cs;^ , and (12)

[0229]

[0230] Di = coeffs; ■ Ci . , „ ..2 1 , (13) wherein A is the covariance vector defined in (10), exp is the exponent matrix (comprising the first row in (9)), coeffs is the coefficient matrix and Di is a vector corresponding to an expectation E[.]. An expectation matrix D is then defined by compounding the D, vectors and, thus, contains in each column an expectation E,[.] - e.g., corresponding to an addend of equation (7) for the two-dimensional example. Matrix D is used to solve Equations (2) and (3) to compute the propagated covariance matrix P,j of the state of uncertainty at time instant L In other words, the uncertainty propagation module 105 determines the covariance and mean of each GME at the time instant f, which is then compounded and output as the GMM of the state uncertainty at time instant

[0231] Afterwards, a conjunction probability module 107 is adapted to compute the conjunction probability rate for the bodies A and B.

[0232] In detail, the conjunction probability module 107 acquires (step 3027):

[0233] a three-dimensional geometry model of the colliding objects, and

[0234] the state uncertainty at time instant t - i.e., {Pi,

[0235] Based on available model of the bodies A and B, the conjunction probability module 107 executes a multi-spheres or a single-sphere solution method (decision step 3029). Advantageously, single sphere solution is chosen when both bodies A and B are modelled as a single sphere, while multi-spheres solution is chosen when at least one of the bodies is modelled by a plurality of spheres.

[0236] If both bodies A and B are modelled as spheres (exit branch S of step 3029), the sphere model of bodies A and B are acquired (step 3031, e.g., from the repository 20) and a plurality of Lebedev quadrature points and weights are computed based on the sphere models (step 3033). In the non-limiting example herein considered, Lebedev's quadrature points are obtained substantially as described in Lebedev, V. (1976): “Quadratures on a sphere" published in USSR Computational Mathematics and Mathematical Physics Vol. 16, No. 2, pp. 10-24.

[0237] If at least one of the bodies A and B is modelled as multi-spheres (exit branch M of step 3029), a spheres mesh modelling of such body A or B and a spheres mesh model or a sphere model of the other body B are obtained (step 3035, e.g., from the repository 20) and plurality of Lebedev quadrature points and weight are computed (step 3037) based on the spheres mesh.

[0238] Particularly, a set of Lebedev points and weights are computed for the multi-spheres model. This set comprises Lebedev points and weights of the spheres forming the external boundary of the multi-spheres model as described below.

[0239] In an alternative embodiment, the repository 20 stores pre-computed Lebedev's points and weights for each one of the bodies monitored. Afterwards, a surface integration is performed to obtain an instant conjunction probability rate (step 3039). In the example considered, in case both bodies A and B are modelled as single spheres surface integration is made over a single hard-body sphere. Conversely, the surface of integration is computed over multiple hard body spheres in case of at least one of the bodies A and B is modelled by a multi-spheres model.

[0240] Particularly, when both bodies A and B are modelled as single spheres, based on the direct method of conjunction probability introduced in Coppola, V. and McAdams, J. V. (2012): “Including Velocity Uncertainty in the Probability of Conjunction between Space Objects" published in A AS12 -247, 22nd, Spaceflight mechanics 2012, Vol. 143, pp.

[0241] 2159-2178, in Akella, M. R. and Alfriend, K. T. (2000): “Probability of Conjunction Between Space Objects" published in Journal of Guidance, Control, and Dynamics Vol. 23, No. 5, pp. 769-772, and in DeMars, K., Cheng, Y., and Jah, M. (2014): “Conjunction Probability with Gaussian Mixture Orbit Uncertainty’ published in Journal of Guidance, Control, and Dynamics Vol. 37, No. 3, pp. 979-985, the conjunction probability rate pc(t) as a function of time is computed as:

[0242] y~ y y i

[0243]

[0244] ■" ■ ";.. .;(14) where:

[0245] R is a joint hard-body radius of bodies A and B (i.e., R = RA + RB),

[0246] Naand Nb are the number of GMEs describing the PDF of bodies A and B respectively,

[0247] Nk is the number of boundary points (in this case equal to the quadrature points) of a surface describing bodies A and B joint together - i.e., a sphere having the radius equal to the sum of the radii of the single spheres approximating bodies A and B, or the shape of the larger between body A and body B, if the smaller body has a negligible size with respect to the larger body,

[0248] wwaflw('>b,t are the weights of the / -th GME of bodies A and B, respectively, at time f,

[0249] ('i)t= - m b,i, with nA, t and m b,i being the mean matrices of the / -th and / -th GMEs of bodies A and B, respectively, at time t,

[0250] A = p<oa f+ with PWa,fand P b,t being the covariance matrices of the / -th and / -th GMEs of bodies A and B, respectively, at time t,

[0251] and v(n) is a relative velocity between bodies A and B defined as:

[0252] — c / f _n_)_eXp J ( - v / .f —n) - ■votii} f, J _erJf ? - v0(n) - 1 L

[0253] I 2;yqn)

[0254]

[0255] 2 [

[0256] where n is a normal vector to a portion dS of the surface S modelling the considered body A or B, while v0(n) and <j2(n) are defined as:

[0257] v0{n) = n7

[0258] and (16)

[0259] <T2fn} = n7' (x^3)1X)'^j n

[0260]

[0261] (17) wherein the joint mean and the covariance matrix are decomposed into position (denoted by the subscript r) and velocity (denoted by the subscript v) components as follows:

[0262] and

[0263]

[0264] Conversely, in case of at least one of the bodies A and B is modelled by a multi-spheres model - e.g., body A, the conjunction probability rate pc(f) is computed as:

[0265] PcfD = 21 y L ^kU'h,uU-Pii UA’WX.yT v(ni.,.tffc)

[0266]

[0267] , (19) where N is a number of boundary points associated with the surface of the of the spacecraft 50 model, i.e. body A, in a radial-transverse-normal (RTN) reference frame referred to the body A, Wbnd.k are the Lebedev weights and Rk are the radii associated with the / r-th boundary point Nk, l^'^1is the relative distance between the body B and the centre of each Lebedev sphere associated with a corresponding boundary point Nk.

[0268] The number of boundary points Nk and Nbnd in equations (14) and (19), respectively, affect the accuracy of the flow integration over the surface of the model of body A or B and, therefore, the overall accuracy of the quadrature method. Generally, the number of quadrature points Nk- which corresponds to the boundary points in the single spheres model, or which comprises the boundary points in the multi-sphere model - is selected based on the relative size of the sphere with respect to the position uncertainty. Moreover, the number of quadrature points Nk linearly affects the computational load of the conjunction probability calculation. The selection of an optimal number of quadrature points Nk depends on specific application requirements and is based on a trade-off between computational load and radius to standard deviation ratio.

[0269] In an embodiment of the invention, a predetermined starting number of quadrature points Nk is selected, preferably based on a contingent scenario considered, and this number is then increased if needed to increase accuracy and / or other constraints e.g., minimum error requirement, computational time available, or due to high bodies complexity. For example, Table 2 reproduced below shows the conjunction probabilities computed for the same conditions but considering a number of quadrature points Nk ranging from 50 to 5810 and shows the relative error between 5810 quadrature points and lower numbers of quadrature points. In this example, the results show that a sufficiently accurate (%Err = 0.0405) computation of the conjunction probability is achieved with approximately as low as 1000 quadrature points.

[0270] Quadrature points 50 590 1000 1454 3074 5810 Collision probability (%} 0.2766552 O.24493O 0.244217 0.2442129 0.2442169 0.2442170

[0271]

[0272] % Error w.r.t 5810 points 13.282 0.291 0.0405 6.5777- 10"42.89- IO"4

[0273] Table 2

[0274] At this point, it is verified if the time instant variable has reached the final instant tf of the conjunction time (decision step 3041). In the negative case (output branch N of step 3041), the time instant variable is increased by a time step At (step 3043) - i.e., a next time instant t+iwill be considered, and the operation is reiterated from step 3007.

[0275] Conversely when the final time instant tf has been reached (output branch Y of step 3041), the conjunction probability Pc for the bodies A and B is obtained by integrating the conjunction probability rate over the conjunction time interval (step 3045). In the considered example, the conjunction probability rate pc(f) is numerically integrated over the conjunction time interval to obtain the total conjunction probability Pc. For example, the numerical integration is performed based on one among a multi-stage, a multi-step or an extrapolation method - e.g., Runge-Kutta (RK), Adams-Bashforth-Moulton (ABM) and Bulirsch-Stoer (BS) respectively. The Applicant found that DOPRI8(7) and Gauss-Legendre (GL-IRK) - belonging to RK algorithms family - and the ABM algorithms provide best results.

[0276] Thus, the conjunction probability module 107 outputs the conjunction probability Pc for the bodies A and B and, preferably, an indication of instantaneous hazard along the surface of one or both the bodies A and B (step 3047). Particularly, the indication of instantaneous hazard comprises computing and outputting a rate of conjunction probability for one or more surface elements of one or more of the bodies considered. The probability conjunction rate is integrated over time and over the entire spacecraft area to obtain the conjunction probability. In other words, this indication provides the knowledge of which part(s) of a body is more susceptible to be hit at each moment in time.

[0277] As anticipated above, the system 1, particularly the conjunction probability arrangement 10, is adapted to pre-compute a respective multi-spheres model for any spacecraft 50, according to a modelling method 5000 of which Figure 7 is a flowchart.

[0278] In general, the modelling method 5000 comprises the following steps.

[0279] Initially, a three-dimensional model of the spacecraft 50 is obtained (step 5001) or, alternatively, computed. Preferably, a simplified three-dimensional model of the spacecraft is computed based on the component division devised in chapter 6 of Chan, F. (2008): “Spacecraft Conjunction Probability' published by Aerospace Press for the method of equivalent cross-section area. For example, a model of the International Space Station model (ISS) consists of 52 elements represented by simple geometrical shapes (parallelepipeds and cylinders) that approximate the ISS at one of its possible configurations, as qualitatively shown in Figure 8.

[0280] Then, a number of spheres Nsph that fit the volume of the spacecraft 50 is determined (step 5003). The number of spheres Nsph depends on the complexity of the shape and the available computational capabilities. Preferably, the spheres radii may change from one sphere to another.

[0281] A sphere-mesh representation is, then, created (step 5005). Particularly, the sphere-mesh representation is defined by a position of the center and a radius of each sphere.

[0282] Preferably, the sphere fitting and mesh generation is computed by means of a sphere mesh method based on spheretree construction described in Bradshaw, Gareth and Carol O'Sullivan, (2004): “Adaptive medial-axis approximation for sphere-tree construction”, published in ACM Trans. Graph., Vol. 23, pages 1-26.

[0283] The sphere mesh method is described below by considering the International Space Station (ISS) as an exemplary spacecraft 50 to be modelled, of which Figure 8 shows a three-dimensional model thereof (selected at previous step 5001).

[0284] The Applicant noted that creating a sphere tree of the solar panels 210 is computationally impractical. Particularly, each panel has a surface area of about 420 m2and a thickness of about 10 cm. By straightforwardly applying the method of Bradshaw and O'Sullivan (2004), more than 50,000 spheres would be required to approximate the parallelepiped shape of each solar panel 210 of the ISS, adding up to a total greater than 420,000 spheres for the set of eight solar panels of the ISS.

[0285] Thus, the sphere mesh method according to an embodiment of the present invention comprises modelling the solar panels 210 and the body 220 of the ISS 200 separately. More generally, the sphere mesh method comprises identifying special portions of the spacecraft 50 - i.e., portions having at least one size substantially smaller (e.g., an order of magnitude lower) than the other sizes. For each special portion a respective sphere-three is computed.

[0286] Particularly, a thicker model of the special portion is created, by increasing the smaller size of the special portion to a desired thickness. In the ISS's example, a thicker model of each solar panel 210 is created, in which each solar panel has a thickness of 10 m. Then, a respective special sphere tree is computed for each thicker model of the special portion {i.e., solar panel 210) and the outer upper and lower boundaries of the expanded sphere tree representing the plane containing the photovoltaic arrays is extracted. The solar panel sphere tree thus obtained is substantially smaller than 50,000 spheres, as qualitatively appreciable in Figure 9, and the complete set of solar panels is described with only 10,000 spheres.

[0287] In series or in parallel, a sphere tree of the remaining portion of the spacecraft 50 is computed. In the ISS's example, the body of the ISS is modelled as a sphere tree using 10,000 sample points and 50,000 spheres by computing a respective body sphere tree. The resulting sphere mesh is illustrated in Figure 10. Despite the complex shape of the body, the sphere mesh produces an accurate replica of the real volume of the ISS. A predetermined number of quadrature points Nk is sampled for each sphere (step 5007). Accordingly, a mesh of Nsph *N points is obtained. A plurality of boundary points Nbnd are extracted from the mesh of Nsph *N points (step 5009). The boundary points Nbnd define an external surface of the volume of the model of the spacecraft 50. As described above, the number of boundary points used is selected based on application-specific constraints and, in general, entails selecting a trade-off between computational time and accuracy.

[0288] Particularly, boundary points of each special sphere tree are extracted and surface points forming the main faces of the special portion are selected. With respect to the solar panel sphere tree, the surface points corresponding to the upper 211 and lower surfaces 212 of the solar panel 210 are selected as qualitatively shown in Figure 9, where the boundary points are bolder than internal points. A boundary extraction algorithm is, preferably, used to identify the external points. In the ISS's example, the upper 211 and lower surfaces 212 are both 5 cm thick slices of the sphere tree. In other words, the upper 211 and lower surfaces 212 represent the quadrature points forming the external boundary of an equivalent parallelepiped with 10 cm thickness. The surface points of the main surfaces are then shifted to the effective size of the special portion ( / .e., the main surfaces are separated by a distance corresponding to the smaller size of the special portion). In the ISS's example, the upper and lower surface points 210 and 220 are shifted to recreate the shape of the solar panel with a thickness of 10 cm.

[0289] In series or in parallel, boundary points of the sphere tree of the remaining portion of the spacecraft 50 are extracted, preferably, by means of a boundary extraction algorithm, as qualitatively shown in Figure 10 where the boundary points are bolder than internal points.

[0290] Once all the boundary points are identified, the boundary points of the one or more special portions (e.g., the solar panels 210) and of the remaining portion (e.g., the body 220) of the spacecraft 50 (e.g., the ISS) are combined as illustrated in Figure 11, obtaining a set of the desired plurality of boundary points NM.

[0291] For each boundary points Nbnd of the plurality of boundary points Nbnd, the associated Lebedev weight Wbnd.k and sphere radius Rk are computed (step 5011) and, preferably, stored in the data repository 20 (step 5013). In other words, Lebedev weights Wbnd,k, sphere centre positions and radii Rk associated with the plurality of boundary points Nbndis the model used to compute the conjunction probability between the spacecraft 50 and another object such as the space debris 60.

[0292] Simulation results

[0293] The accuracy of a conjunction probability calculation method according to the present invention is assessed by comparing the results with Monte Carlo simulations. To this extent, a database of satellite conjunction with Monte Carlo analysis reported in Alfano, S. (2009): “Satellite conjunction Monte Carlo analysis" published in Advances in the Astronautical Sciences Vol. 134, pp. 2007-2024, is used for verification. Alfano studied 12 different cases for close conjunctions at LEO, MEO, HEO and GEO with 108 Monte Carlo samples and considering only the central gravitational acceleration.

[0294] Particularly, cases 7 and 12 present special properties and are therefore selected. Case 7 involves the relative motion between two satellites in LEO which results in a conjunction probability of Pc ~ 10-4. This allows to evaluate and compare the accuracy of the method to predict conjunctions with a very low risk. It must be noted that 109 Monte Carlo samples are required to compute this probability with sufficient accuracy. Case 12 presents a peculiar scenario of two satellites co-located in identical LEO orbits and with the same initial state uncertainty. The propagated state and uncertainty will be identical for both satellites, but a conjunction does not necessarily occur. This is a special event since the conjunction develops continuously.

[0295] Table 3, below, presents the initial position and velocity of the satellites involved in the two cases 7 and 12 selected.

[0296] Initial Position x [km] y [km] z |km]

[0297] Satellite A 8397.9865385122 1889.3098033002 1889.3098033002

[0298] Case 7

[0299] Satellite B 6938.8785812399 1888.2177447937 1888.1503815517

[0300] Case 12 Satellites A & B 8878.137 8.0 M

[0301] Initial Velocity S [km] y [km] i ]km]

[0302] Saieltfte A -2.9571994197397 4.9601786006789 4.9601786006789

[0303] Case 7

[0304] Satellite B -2.9552934327179 4.9608101452466 4.9606345630827

[0305]

[0306] Case 12 Satelli tes A & B 0,0 7.6126081732239 0.0

[0307] Table 3

[0308] The initial state uncertainty is presented in following Table 4. It must be noted that the covariance elements excluded from the table are set to zero and that the covariance matrix is symmetric ( / . e., ox,y= oyx). Covariance Case 7, Sat A Case 7, Sat B €ase 12, Sat A «t B ttiits <!.409t!7iS84t!37ffS ■ nr43.4IO728283379I • It!"48.57B2747<X«mS9<f - Iff"332.3460130768548- IO"'12.3335737435007 ■ IS"44.421725299O5SI ■ 1O-S2.34i)0iff«7B&i4S- -O"42.3386979731152 - Iff"44.8- :O"i;

[0309] i fair s 0.8960443775773 ■ -O”39.8916238832886- 10"" -1 3395230358788 IO"'3S.89G844377577.3- 1O"59.8912538411159 • it!"40.6

[0310] -1.8508869231852 - IB"4-5 .6603641405277- IO"48.0

[0311] 8.5660758872296 10"u3.5079630985814- 10”i s5.8- 10"14

[0312] 5.79S&6495S8&52 • W115.7958556778401 ■ 10”u1.8- it!"34

[0313] 5.7989640563652 ■ Iff”3 15.791662234785-; :i.O- Iff "34S km- $ 2.5<>57S®Sa3aBM • IO"4* 2.55145288093237- 16"n0.6

[0314] 2.5057999830694 - 10"112.5044351557813 ■ 1ft”3 80.0

[0315]

[0316] -4.2038359436148 - Iff”’1-4.2B3991946S489 • 16"118.0

[0317] Table 4

[0318] In the simulations run by Alfano (2009), the conjunction time is set from TCA -1420 s to TCA+1420 s. To be consistent with the results from the Monte Carlo analysis, the same time span is selected for the conjunction probability calculation performed by the method of the present invention. By considering that the conjunction probability must be integrated over an extended period, the simulations are performed with a Taylor expansion order k of three (k = 3). Further, the main drivers for the computational load on the Pc calculation process using the time integration method are the number of GMEs and the number of integration steps (which is a function of the conjunction period and of the time step). In a first evaluation a maximum of 51 GMEs is selected for the GMM. In a second test, a maximum of 37 GMEs is performed to evaluate the sensitivity of the method. Finally, the time step is set to 10 seconds.

[0319] It should be noted that the method according to the present invention not only provides the total conjunction probability, but also the evolution of the conjunction probability rate and the cumulative conjunction probability over time. This provides substantial additional information than the risk as a single number as provided by the method known in the art.

[0320] Figures 12A and 12B show the evolution of the conjunction probability rate and the conjunction probability for Case 7 with a lead time of 48 hours is presented. As readily observed, the real duration of the conjunction is only about 300 seconds, which is relatively small with respect to the 2840 second window studied. Figures 13A and 13B show the evolution of the conjunction probability rate and the conjunction probability for Case 12. In this case, the probability rate of conjunction spans through the entire time studied, which is expected since both satellites follow a nominal trajectory. Although the conjunction probability obtained for Case 12 is one order of magnitude higher than that of Case 7, the probability rate is in the same order of magnitude in both cases. This hints that the total conjunction probability alone - as provided by the systems at the state of the art - might be insufficient to correctly determine conjunction risks in some scenarios. The results obtained are verified with the results from a Monte Carlo analysis, see Alfano (2009, Fig A7, Fig A12).

[0321] Table 5 presents the resultant conjunction probability calculated with the method of the present invention (indicated as DA-GMM) and compared to alternative methods provided by Alfano (2009). In the table the cases where the method provides and improved computed conjunction probability than the alternative methods are highlighted. For case 7, an error of 0.04% is achieved, which shows an improvement with respect to both the linear and non-linear alternative known methods. Method PrAPei%) PcAPC(%)

[0322] Case 7 Casel 2.

[0323] Monte Carlo samples) 1.51 18"4-6.48 2,55 ■ 18-38.23 Maute Carks {ltT!samples} KOI 462 - IO"42,55595 - IO"3

[0324] DA-GMM (5.1 GMEs) 1.615292 - 1©48.84 2.444178 - IO"34.37 DA-GMM(37GMEs) LW&687- IO'"'5-0.34 2.443430 -8 ’ 4,48 Voxels » = 50 1.64414- 18-41.33 3.636683- SO"343.07 Voxels n = 180 1.617 IS -IO”40.16 3.079178- 18“’3

[0325] Adjoining cylinders » = 58 1.61677 - IO"40.13 8.9 100. ©0 Adfoloing tenders n. = 10© 1.61677- IO”40.13 0.8 100.08 Parallelepipeds e = 50 1.63761 - IO"40.18 8.0 100.00 Parallelepipeds n = 100 1.61701 -10“40.15 0.8 180.00 Linear Alfano n = 58 L58147- 18"4- 2.05 1.817488- I0"s-24.98 Linear Alfeno u = 108 1.58147 -IO"4-2.05 1.917487- 18"’3-24.88 Linear Patera H = 50 L58146-18"4-2.05 1.817487- 10"s-24.98

[0326]

[0327] Linear Patera n = 180 1.58146- IQ"4-2.05 1.917487- 18"’3-24.88

[0328] Table 5

[0329] Moreover, the total time to compute the conjunction probability over the 2840 second interval with a time-step of 10 seconds is three hours and less than one hour for the cases with 51 and 37 GMEs, respectively. This allows to compute the conjunction probability with sufficient time to conduct an avoidance manoeuvre and therefore it can be used during operations.

[0330] For case 12, the results obtained by both using 37 and 51 GMEs show a significant improvement over the known methods. The inferior performance of the known methods is due to the fact that there is no relative velocity, leading to a zero cumulative probability. A similar effect is produced in scenarios with low relative velocities, such in formation flying. For these scenarios, the known non-linear methods fail, and the known linear methods provide a low accuracy due to the large conjunction times. Moreover, these known methods under-estimate the conjunction probability, which can lead to unexpected casualties. On the contrary, the method according to the present invention estimates the conjunction probability with an error of 4.37%, for Case 12, providing the closest approximation to the real magnitude of the risk w.r.t. the other known methods.

[0331] In summary, the method according to the present invention features the following advantages over the prior art: • a conjunction probability error of 0.04% is achieved with a lead time of 48 hours,

[0332] • the method is not limited to cases with large relative velocity,

[0333] • a high accuracy estimate can be obtained with 37 GMEs,

[0334] • low conjunction probabilities (Pc ~ 10-4) are correctly predicted,

[0335] • the accuracy in the conjunction probability Pc calculation is improved by a factor greater than 70% with respect to the best-known method, and

[0336] •an accurate conjunction probability computation for the three-dimensional geometry, applicable to any conjunction geometry, is obtained within five minutes for one GME and within one hour for 51 GMEs, which makes the method particularly apt for operational conjunction avoidance.

[0337] The invention thus conceived is susceptible to numerous modifications and variations, all of which are within the scope of the inventive concept that characterizes it.

[0338] For example, in a different embodiment of the system, the monitoring assembly and / or the mission control assembly directly provides data to the conjunction prediction assembly. In this case, the conjunction prediction assembly comprises at least one memory element adapted to store such data.

[0339] As an addition or an alternative, the multi-sphere models are pre-computed by data processing assembly - for example, a computer network - and then stored in the repository or in the conjunction prediction assembly.

[0340] In an embodiment, the conjunction prediction assembly is managed independently from the mission control assembly and / or the monitoring assembly. In such a case, conjunction predictions for one or more spacecraft can be requested on demand by one or more respective mission control assemblies by establishing a communication channel with the conjunction prediction assembly through which messages are exchanged; such messages containing, in a non-limiting manner, conjunction prediction request, conjunction predictions, spacecraft-related data and / or conjunction data messages.

[0341] Nothing prevents modeling bodies other than spacecrafts, such as for example a debris or a celestial body via the multi-spheres modelling method of above mutatis mutandis.

[0342] In an embodiment of the invention, the time step used to perform the computation is defined as follows. The time step is generally a case-dependent variable, which is selected as a function of the shape of the probability rate curve and the integration time. In case of very short term conjunctions, a shorter time step will be required. For short term conjunctions (typically conjunction time interval of less than five seconds) a time step equal to 0.1 second is preferably used as recommended by Coppola, V. and McAdams, J. V. (2012): “Including Velocity Uncertainty in the Probability of Conjunction between Space Objects" published in AAS12247, 22nd, Spaceflight mechanics 2012, Vol. 143, pp. 2159— 2178.

[0343] On the contrary, in case of long term conjunctions the conjunction probability is spread over a longer conjunction time interval, which allows selecting longer time steps, preferably, in the order of 10-60 seconds.

[0344] The Applicant has found that the time step selected to have 100-300 integration points - i.e., time instants, also indicated as time steps, in which is performed the integration of formula (15) - within the conjunction time interval provides optimal results. Optionally, the time step is varied based on the velocity of the bodies in fast velocity changing trajectory segments if required.

[0345] In general, the number of time steps depends on the ratio between the conjunction time interval and the time step. Nonetheless, in an embodiment of the invention, a minimum number of time step (e.g., a highly conservative number such as 10000 time steps) can be arbitrarily selected regardless of the length of the conjunction time interval.

[0346] In an embodiment of the invention, the hard-ball radius - i.e., the radius of the sphere to which a body A or B is modelled - depends on the geometrical specifications of the problem, particularly, the size of the body A or B. However, due to the approximation required to calculate the surface integration along the sphere with Lebedev's method, the method will fail when the sphere is too big compared to the position uncertainty. Thus, the hard-ball radius of the body A or B is selected at least one order of magnitude smaller than the standard deviation of the position uncertainty of the body's A or B.

[0347] Generally, the DA-GMM model is based on the assumption that by dividing the initially Gaussian distribution into a large set of sub-distributions which are also Gaussian, the propagation of the individual elements through the nonlinear dynamics will correctly approximate the final non-Gaussian distribution. This implies that the error in keeping the individual GMEs as Gaussian distributions through the propagation must be negligible. Accordingly, in an embodiment of the invention, the conjunction prediction assembly comprises a validation process to verify that the cumulative error associated with the GMEs of the DA-GMM model remains sufficiently low - e.g., lower than an error threshold such as 10'8. Furthermore, this validation process is used also to select the number of required GMEs.

[0348] Preferably, the validation process comprises computing an L2 error between approximations having different GMEs to evaluate the effect of increasing GMEs on the final distribution. Then the validation process comprises a goodness of fit test procedure, which is adapted to test if a given data set is normally distributed. Non-limiting examples of suitable goodness of fit tests comprise Pearson's x2test, the Kolmogorov-Smirnov (KS) test and the Anderson-Darling (AD) test. The information obtained through the validation process allows defining a trade-off between computational time and accuracy.

[0349] In an embodiment of the invention, once the conjunction probability rate is computed, it is checked whether a probability rate threshold is violated and, in the affirmative case, the conjunction prediction assembly performs one or more additional simulations leveraging model parameters in order to comply with the probability rate threshold.

[0350] As shall be clear to the skilled person, one or more method steps or system elements can be replaced by other technically equivalent ones, falling within the scope of the attached claims, according to the specific requirements. For example, nothing prevents that, in a simpler embodiment of the method, the conjunction probability computation based on DA-GMM is based only on single sphere models of the bodies. Conversely, nothing prevents that in different embodiments the multi-sphere models obtained as described above are used to compute a conjunction probability according to a different method (e.g., known methods that output a propagated covariance - such as state transition tensor (Romgens et al., 2011), polynomial chaos etc.).

[0351] As shall be apparent to those skilled in the art, the methods and systems according to the present invention can be used to compute a conjunction probability between two spacecrafts (instead of a spacecraft and space debris), as well as any other pair bodies in space, mutatis mutandis.

[0352] Furthermore, it should be noted that one of the bodies considered for assessing conjunction probability, may consist of a cluster of objects, e.g. a cluster of space debris having substantially corresponding relative velocities w.r.t. the other body considered and being within a predetermined maximum distance one from the others.

[0353] In a different embodiment of the invention, the system 1 is adapted to monitor and perform a conjunction-reaching manoeuvre, also indicated as a rendez-vous manoeuvre. In other words, the systems are adapted to monitor and regulate proximity manoeuvres during a rendez-vous between two bodies in space - e.g., two spacecraft - to maximise the probability of conjunction - i.e., the conjunction probability Pc is brought as close as possible to one. In this case, the system 1 performs a method for rendez-vous 6000, of which Figure 14 is a flowchart. Method 6000 comprises steps 6001 to 6007 that are substantially equivalent to steps 1001 to 1007 of method 1000 described above and are herein not repeated for the sake of brevity.

[0354] After step 6007, the conjunction prediction assembly 10 determines if one or more of the computed conjunction probabilities is below a low conjunction threshold (decision step 6009). For example, the low conjunction threshold is equal to a conjunction probability equal to or lower than 10%.

[0355] If there is not any conjunction probability below the low conjunction threshold (output branch N of step 6009), the procedure is interrupted, since there is not any risk of missing the rendez-vous. Accordingly, the conjunction prediction assembly 10 enters an idle state, waiting for a data update and / or a start command to reiterate the method from step 6001.

[0356] If one or more of the conjunction probabilities are below the low conjunction threshold (output branch Y of step 6009), for each conjunction probability below the conjunction threshold a corresponding conjunction data message is generated (step 6011). In this case the data, the generic conjunction data message comprises information identifying and / or defining the spacecrafts 50 associated with the warning, telemetry information - e.g, Time of closest approach (TCA), conjunction data message (CDM) creation date, trajectory, relative speed, etc. of the spacecrafts 50.

[0357] The conjunction data messages are provided to the mission control assembly 30 (step 6013), which in its turn analyses each conjunction data message to determine if a rendez-vous adjustment manoeuvre is necessary (decision step 6015). In the negative case (output branch N of step 6015) there is not (yet) any need for a rendez-vous adjustment manoeuvre and the method comprises waiting for a data update and / or a start command to reiterate the method from step 6001.

[0358] On the contrary, if a rendez-vous adjustment manoeuvre is required (output branch Y of step 6015), the mission control assembly 30 determines the parameters characterizing the rendez-vous adjustment manoeuvre (step 6017) to be performed by at least one of the spacecrafts 50. Then, the mission control assembly 30 send manoeuvre instructions to the selected spacecraft(s) 50 to perform the rendez-vous adjustment manoeuvre (step 6019) and monitors the execution thereof to its completion (step 6021). For example, the mission control assembly 30 stores one or more instructions sets that allows implementing one or more procedures to assess the necessity of rendez-vous adjustment manoeuvre, to compute the rendez-vous adjustment manoeuvre based on the conjunction data messages and to command the spacecraft 50 to execute the rendez-vous adjustment manoeuvre.

[0359] Afterwards, the method 6000 proceeds to wait for a data update and / or a start command to reiterate the method from step 6001.

[0360] The improvements provided by the present invention are not limited to the aerospace industry. In different embodiments, the method according to the present invention is adapted to predict conjunctions between bodies in Earth's atmosphere, on land or water such as aircrafts, ground vehicles or naval vessels both manned and unmanned. An embodiment of the present invention schematically shown in Figure 15, comprises a conjunction prediction system 1 A adapted to implement the method of conjunction prediction to estimate the probability of collision between aircrafts 50A, and plan trajectory adjustments based upon the probability of collision. For example, the system is implemented to improve safety in regions that are not covered by any active Air Traffic Control (ATC) infrastructure, but controlled entry to such regions by ATC infrastructure exist. This is the case of uninhabited open spaces, such as the oceans, and restricted airspaces (e.g., Afghanistan's airspace).

[0361] In another embodiment, an embedded conjunction prediction system 1B comprises an onboard conjunction prediction assembly 10B, which cooperates with an autonomous or semi-autonomous onboard flight management system 51 of at least one of the aircrafts 50B, as schematically shown in Figure 16. Preferably, the conjunction prediction assembly 10B maintains a communication channel - e.g., via radiofrequency - with one or more among the data repository 20 B, the mission control assembly 30 and the monitoring assembly 40. Even more preferably, the conjunction prediction assembly 10B receives data, e.g. flight data, from the flight management system 51 or from another onboard control system (known in the art and not shown for the sake of simplicity).

[0362] The system 1A, 1B executes the methods 1000A-6000A of Figures 17-22, corresponding to methods 1000-6000 described above mutatis mutandis. In particular, the following differences are implemented. In these embodiments, the collision avoidance manoeuvre is provided to avoid an object, e.g. another aircraft, and the rendez-vous manoeuvre is provided to approach a body, e.g. to perform operations such as in-flight refueling, intercept mission, etc.

[0363] The data stored in the data repository 20 comprises data relating to aircrafts, e.g. 3D models of monitored aircrafts. The DA module 103 and the propagation module 105 are adapted to compute the dynamics of aircrafts, rather than spacecrafts. For example, the DA module 103 receives aerodynamic model inputs, reference frames, and body parameters consistent with conventional 6-DOF aircraft dynamics according to standard aircraft flight dynamics models as known in the art., rather than gravitational perturbations. Further, the propagation module 105 is configured to handle aerodynamic flight dynamics consistent with standard aircraft navigation and estimation techniques known in the art -instead of orbital motion.

[0364] In addition, while the GMM module 101 remains substantially unchanged, an appropriate source for the last known states and corresponding last state covariance of the aircrafts is connected to the system (to perform step 2001A). These input data can be computed and provided based on ATC data, GNSS data. For example, ground-based tracking systems, such as Primary and Secondary Surveillance Radar (PSR and SSR),— comprised in the monitoring assembly 40 - and / or Automatic Dependent Surveillance-Broadcast (ADS-B) provide position and velocity measurements of aircrafts provide data that are utilisable to compute state and associated covariance of aircrafts 50A, 50B in their operating range. Additionally, Unscented Kalman Filters (UKF) and / or Extended Kalman Filters (EKF), typically executed by the onboard flight management system 51 are usually fed with a combination of GNSS data, Inertial Navigation System (INS) data, sensors measurements to output the aircraft 50B state and associated covariance. In case of a manned aircraft, optionally, the outputs provided by the methods 1000A and 6000A comprise manoeuvre instruction that are provided to a pilot of at least one of the aircrafts trough an aircraft interface for performing the collision avoidance manoeuvre or the rendez-vous manoeuvre, as an alternative or as an addition to the command for performing the same manoeuvres provided to the onboard control system of the aircraft(s).

[0365] As schematically shown in Figures 23 and in Figures 31, a conjunction prediction systems 1C, 1E is adapted to predict a conjunction between two bodies on land, e.g. two ground vehicles 50C, or two bodies in water, e.g. two naval vessels 50E, respectively. Similarly, conjunction prediction systems 1D, 1 F comprise an onboard conjunction prediction module 10D, 10F embedded in at least one ground vehicle 50D or at least one naval vessel 50F, respectively, as schematically shown in Figures 24 and 32. The onboard conjunction prediction modules 10D, 10F are adapted to exchange data with an autonomous or semi-autonomous onboard piloting system 51 of the respective ground vehicle 50D or naval vessel 50F. The systems 1C-1 F have features corresponding to the features of system 1 described above unless otherwise noted.

[0366] The systems 1C-1D and 1E-1F execute methods 1000C-6000C and 1000E-6000E, of which Figures 25-30 and 33-38 are flowcharts, respectively. These methods correspond to methods 1000-6000 described above mutatis mutandis. In particular, the following differences are implemented. In these embodiments, the collision avoidance manoeuvre is provided to avoid a body, e.g. another ground vehicle or naval vessel, a natural or artificial structure etc., and the rendez-vous manoeuvre is provided to approach a body, e.g. to perform operations such as folding in line, maintaining a formation docking, contacting, etc.).

[0367] The data stored in the data repository 20 comprises data relating ground vehicles or naval vessels according to the embodiment considered, e.g. 3D models of ground vehicles or naval vessels. The DA module 103 and the propagation module 105 are adapted to compute the dynamics of ground vehicles 50C,D or naval vessels 50E,F, rather than spacecrafts. For example, the DA module 103 receives aerodynamic model inputs, reference frames, and body parameters, full rigid-body equations, hydrodynamic damping, environmental effects, etc. consistent with conventional 6-DOF dynamics for ground vehicles or naval vessels according to standard ground vehicles or naval vessels travel dynamics models as known in the art. rather than gravitational perturbations. Further, the propagation module 105 is configured to handle aerodynamic flight dynamics consistent with standard ground vehicles or naval vessels navigation and estimation techniques known in the art.

[0368] In addition, while the GMM module 101 remains substantially unchanged with respect to the previous embodiments, an appropriate source for the state covariance of the ground vehicles or naval vessel is connected to the system 1C-1F (to perform step 2001C-D). For example, ground-based tracking systems, such as radar systems, plate reader system etc. or on board system such as GNSS, automatic identification system (AIS), sensors that provide position and velocity measurements can be used to compute state and associated covariance of at least one of the ground vehicles 50D or naval vessel 50F - e.g., the computation of the state and associated covariance might be performed by a control unit 51 of the ground vehicle 50D or of the naval vessel 50F or by an additional module (not shown) embedded in the ground vehicle 50D or the naval vessel 50F.

[0369] In case of a manned ground vehicle or naval vessel, optionally, the outputs provided by the methods 1000C.D and 6000C.D comprise manoeuvre instruction that are provided to a pilot of at least one of the ground vehicles or naval vessels trough an interface thereof for performing the collision avoidance manoeuvre or the rendez-vous manoeuvre, as an alternative or as an addition to the command for performing the same manoeuvres given to the onboard control system of the ground vehicle(s) or naval vessel(s).

[0370] In further embodiments of the invention, the methods and systems of above are configured to cooperate with a defense system to determine collision trajectories between ammunition and a moving target such as a missile or an unmanned aircraft.

[0371] In an exemplary embodiment, an impact prediction system 1G, 1H is connected to a control system 52 of an artillery piece 50G or is embedded in a ballistic missile 50H, as schematically shown in Figures 39 and 40.

[0372] The data stored in the data repository 20 comprise data relating intercepting means 53, 50H, and target T e.g. a 3D models of the relevant intercepting means 53, 50H and target T. The DA modules 103 is fed input data and the propagation module 105 is modified in such a way that the dynamics reflect the equations of motion of the assessed intercepting means 53, 50H, and the environment in which the projectile travels - i.e.: intercepting means starting speed / acceleration data, propulsions data, atmosphere, space or water dynamics parameters in a manner similar to what described with respect to systems for spacecrafts, aircrafts, ground vehicles and naval vessel mutatis mutandis. Preferably, the missile 50H comprises an onboard conjunction prediction assembly 10H connected to a mission control module 51 of the missile 50H, which operates in a similar manner as described above with respect to onboard conjunction prediction assembly 10B, 10D and 10F described above. Alternatively, the conjunction prediction assembly 10 on the ground is used.

[0373] In case of a moving target T, data relating to the intercepting means such as a shell 53 or the missile 50H, and of the target T are acquired, and the impact prediction system 1G, 1H is configured to perform the most suitable set of methods among methods 2000-5000, 2000A-5000A, 2000E-5000E, mutatis mutandis, according to the environment in which the intercepting means travels.

[0374] Further, in the case of the missile 50H one among the methods 1000, 6000, 1000A, 6000A, 1000E, 6000E is performed, for adjusting the trajectory of the missile 50H to avoid hitting a specific body (methods 1000, 1000A, 1000E) or to increase the probability of hitting the target T (methods 1000, 1000A, 1000E).

[0375] Conversely, in case of the artillery piece 50G a method 7000, of which figure 41 is a flowchart, which differs from method 6000 in what follows. For each conjunction data message informing the control assembly 30 that the conjunction probability is below the low conjunction threshold (decision step 7015), the control assembly determines a gun laying adjustment (step 7017) and commands the artillery piece to perform such adjustment (step 7019) and, monitors the adjustment (step 7021). Conversely, if the CEP is equal or lower than the maximum CEP threshold (output branch N of step 7015), the control assembly 30 commands firing the shell (step 7023).

[0376] In case of a static target T, the conjunction prediction and target reaching system 1G, 1H execute methods 2000G-5000G, Figures 41-46, corresponding to methods 2000-6000 described above mutatis mutandis. In particular, the following differences are implemented.

[0377] Particularly, data relating to the intercepting means such as an artillery shell 53 or the missile 50H, and of the target T are acquired. Since the target is static, methods steps 2003G - 2007G of method 2000G comprise computing GMEs and GMM at the initial state (i.e., launch) only of the intercepting means 53, 50H are evaluated. Further, in method 3000H steps 3001H-3039H are performed in such a way to compute the final covariance matrix of the impact position of the intercepting means only (at step 3039G). Afterwards, (at step 3045G), a Circular Error Probable (CEP) is computed based on the final impact covariance matrix. For example, the CEP is expressed as a function of the standard deviations of the positional uncertainties along the principal axes based on Nelson, W. (1988): "Use of Circular Error Probability in Target Detection”, The MITRE Corporation, ESD-TR-88-109, Bedford, MA.

[0378] Again, in the case of the missile 50H one among methods 1000G and 6000G substantially corresponding to methods 1000, 6000, mutatis mutandis, is performed, for adjusting the trajectory of the missile 50H to avoid hitting a specific body (method 1000G) or to increase the probability of hitting the target T (method 6000H).

[0379] Conversely, in case of the artillery piece 50G a method 7000G, of which Figure 48 is a flowchart, which differs from method 7000 in what follows. For each conjunction data message informing the control assembly 30 that the CEP is higher than a maximum CEP threshold (at decision step 7015G), the control assembly determines a gun laying adjustment (step 7017G) and commands the artillery piece to perform such adjustment (step 7019G) and, monitors the adjustment (step 7021 G). Conversely, if the CEP is equal or lower than the maximum CEP threshold (output branch N of step 7015G), the control assembly 30 commands firing the shell (step 7023G).

Claims

CLAIMS1. Method (1000-6000; 1000A-E - 6000A-E) of predicting a conjunction between two bodies (50, 60; 50A,B; 50C,D; 50E,F, 50H) comprising the steps of:acquiring (2001; 2001 A-G - 2001 A-G) data defining at least one last known state and last known covariance of the two bodies;computing (2003-2007; 2003A-G - 2007A-G) a plurality of Gaussian Mixture Elements (GMEs) based on the acquired data;computing (3001-3003; 3001A-G - 3003A-G) a Gaussian Mixture Model (GMM) of an initial state based on the GMEs, the initial state being a state of the bodies at an initial time instant defining a lower boundary of a conjunction time interval comprising a time of closest approach of the two bodies;computing (3005-3015; 3005A-G - 3015A-G) Taylor expansions of dynamic equations of motions of the bodies for each GME in the GMM;computing (3017-3025; 3017A-G - 3025A-G) a final state and a GMM of the final state of the two bodies based on the Taylor expansions, the final state being a state of the bodies at a final time instant defining a higher boundary of the conjunction time interval;computing (3039; 3039A-G) a conjunction probability rate as a function of time based on the initial and final states and GMMs of the initial and final states, andcomputing (3045; 3045A-G) a conjunction probability by integrating the conjunction probability rate over the conjunction time interval, anddetermining (1015; 1015A-G) a possible conjunction if the conjunction probability exceeds a predetermined threshold, characterized in thatwherein the step of computing (3017-3025; 3017A-G - 3025A-G) a final state and a GMM of the final state of the two bodies comprises:defining (3019; 3019A-G) a first matrix and a second matrix based on the Taylor expansion of state dynamic equations, the first matrix comprising coefficients of the Taylor expansion and the second matrix comprising exponents of the Taylor expansion,defining (3021; 3021 A-G) a covariance matrix based on the GMM of the initial state, andcomputing (3023; 3023A-G) the GMM of the final state based on a combination of the first matrix, the second matrix and the covariance matrix.

2. Method (1000-6000; 1000A-G - 6000A-G) according to claim 1, wherein the step of defining (3021; 3021A-G) a covariance matrix comprises:defining a vector of initial state covariance values comprised in the GMM, anddefining the covariance matrix by compounding a plurality of vectors of the initial state covariance values, preferably the number of vectors comprised in the covariance matrix being equal to the maximum number of permutations of initial state covariance values.

3. Method (1000-6000; 1000A-G - 6000A-G) according to claim 1, wherein the step of computing (3023; 3023A-G) the GMM of the final state based on a matrix combination comprises computing an expectation operator matrix as:D — CO6ffSxC[ixn^comp]where D is the expectation operator matrix, coeffs is the first matrix, ncompis the number of state variables, a C is a matrix defined as:where npermsis the maximum number of permutations of initial state covariance values, and B is a matrix defined as:B =4exp11x / nslwhere A is the covariance matrix and exp is the second matrix.

4. Method (1000-6000; 1000A-G - 6000A-G) according to claim 3, wherein the step of computing (3023; 3023A-G) the GMM of the final state based on a matrix combination comprises:computing a covariance matrix P / j of the final state as:Pi) °-V.’ j IW)where Ci,pi...pnand Cj,qi...qnare coefficients of the Taylor expansion, IJIJ is the mean of the / -th and / -th components of a final state Xf, and 5xf1+(<1and 6xr”*cinare deviations, and E is the expectation operator, andcomputing the mean , of the / -th component of a final state as:J f?s +~ J5. Method (1000-6000; 1000A-G - 6000A-G) according to any one of the preceding claims, wherein the step of computing (2003-2007; 2003A-G - 2007A-G) a plurality of GMEs comprises:selecting (2003; 2003A-G) a number of GMEs to be computed;selecting (2005; 2005A-G) a univariate splitting library based on the number of GMEs, andscaling (2007; 2007A-G) the univariate splitting library to fit a multivariate Gaussian distribution, defined by the initial states and by the initial state covariance of the two bodies.

6. Method (1000-6000; 1000A-G - 6000A-G) according to any one of the preceding claims, wherein the step of computing (3023; 3023A-G) a final state and a GMM of the final state of the two bodies comprises:computing a Taylor expansion of the final state for a deviation from the initial state for each GME of the GMM, and computing a Taylor expansion of the final state for a deviation from the initial state for the GMM by combining the Taylor expansions computed for each GME.

7. Method (1000-6000; 1000A-G - 6000A-G) according to any one of the preceding claims, wherein the step of computing (3005-3015; 3005A-G - 3015A-G) the Taylor expansions of dynamic equations of motions of the bodies for each GME in the GMM comprises exploiting a Differential Algebra Computational Engine algorithm, which for each body is configured to:defining the initial state of the body as a differential algebra variable in an inertial reference frame;calculating environment variables as a function of the initial state;defining aerodynamic forces, associated to the environment variables, in the inertial reference frame; computing an aerodynamic acceleration of the body based on the aerodynamic forces and properties of the body; computing gravitational acceleration due to Earth's gravity by determining respective Legendre polynomials and computing potential gradient in spherical coordinates;computing a total acceleration by adding the previously-computed accelerations, andnumerically integrating, by means of an integrator based on Runge-Kutta-Fehlberg (RKF).

8. Method (1000-6000; 1000A-G - 6000A-G) according to any one of the preceding claims, wherein the step of computing (3005-3015) the Taylor expansions of dynamic equations of motions of the bodies for each GME in the GMM comprises exploiting a Differential Algebra Computational Engine algorithm, which for each body is configured to:defining the initial state of the body as a differential algebra variable in an inertial reference frame;calculating environment variables as a function of the initial state;defining aerodynamic forces, associated to the environment variables, in the inertial reference frame; computing an aerodynamic acceleration of the body based on the aerodynamic forces and properties of the body;computing gravitational acceleration due to Earth's gravity by determining respective Legendre polynomials and computing potential gradient in spherical coordinates;computing third body acceleration from Sun and Moon based on the positions of Sun and Moon retrieved from ephemeris in a SPICE library;computing a total acceleration by adding the previously-computed accelerations, andnumerically integrating, by means of an integrator based on Runge-Kutta-Fehlberg (RKF).

9. Method (1000-6000; 1000A-G - 6000A-G) according to any one of the preceding claims, wherein the step of computing a conjunction probability rate comprises selecting between a single sphere model and a multi-spheres model of each body, andperforming a surface integration based on the selected model for each body,wherein at least one multi-sphere model is computed by:acquiring a three-dimensional model of the body,fitting a plurality of spheres inside the three-dimensional body,generating a sphere-mesh representation based on the plurality of spheres,extracting a plurality of boundary points from the sphere-mesh representation, the boundary points defining a surface of the body, andcomputing Lebedev weights, sphere centre positions and radii associated with the plurality of boundary points selected.

10. Method (1000-6000; 1000A-G - 6000A-G) of preventing a collision between at least two bodies in space, the method comprising:executing the method of predicting a conjunction between two bodies in space according to any one of the preceding claims, andif a possible conjunction is determined by the method of predicting a conjunction:generating a corresponding conjunction data message;computing a conjunction avoidance manoeuvre based on the conjunction data message, and command at least one of the two bodies in space to perform the conjunction avoidance manoeuvre.

11. Method (1000-6000; 1000A-G - 6000A-G) of performing a rendez-vous between at least two bodies in space, the method comprising:executing the method of predicting a conjunction between at least two bodies in space according to any one of the preceding claims 1-9, andif a possible conjunction is determined by the method of predicting a conjunction:generating a corresponding conjunction data message;computing a rendez-vous manoeuvre based on the conjunction data message, andcommand at least one of the two bodies in space to perform the rendez-vous manoeuvre.

12. System (1) of automated conjunction prevention between at least two bodies in space comprising:a conjunction prediction assembly (10), comprising at least a computing element, which is adapted to predict a conjunction probability between at least two bodies in space by implementing the method according to claim 10, and a mission control assembly (30) which is adapted to compute, based on the conjunction prediction, and command to perform a conjunction avoidance manoeuvre to at least one body of the two bodies in space.

13. System (1) of automated rendez-vous between least two bodies in space comprising:a conjunction prediction assembly (10), comprising at least a computing element, which is adapted to predict a conjunction probability between at least two bodies in space by implementing the method according to claim 11, and a mission control assembly (30) which is adapted to compute, based on the conjunction prediction, and command to perform a rendez-vous manoeuvre to at least one body of the two bodies in space.

14. System (1A,1B) of conjunction prediction and reaction for at least one aircraft (50A; 50B) comprising: a conjunction prediction assembly (10; 10B), comprising at least a computing element, which is adapted to predict a conjunction probability between the aircrafts and a further body within earth atmosphere by implementing the method according to claim 10 or 11, anda control assembly (30) which is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the aircraft a conjunction avoidance manoeuvre to be performed or to command a control system of the aircraft to perform the conjunction avoidance manoeuvre, orwherein the control assembly (30) is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the aircraft a rendez-vous manoeuvre or to command a control system of the aircraft to perform the rendezvous manoeuvre.

15. System (1 C; 1 D) of conjunction prediction and reaction for at least one ground vehicle (50C; 50D) comprising: a conjunction prediction assembly (10; 10D), comprising at least a computing element, which is adapted to predict a conjunction probability between the ground vehicle and a further body on land by implementing the method according to claim 10 or 11, anda control assembly (30) which is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the ground vehicle a conjunction avoidance manoeuvre to be performed or to command a control system of the ground vehicle to perform the conjunction avoidance manoeuvre, orwherein the control assembly (30) is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the ground vehicles a rendez-vous manoeuvre or to command a control system of the ground vehicle to perform the rendez-vous manoeuvre.

16. System (1 E; 1 F) of conjunction reaction between for at least one naval vessel (50E; 50F) comprising: a conjunction prediction assembly (10; 10F), comprising at least a computing element, which is adapted to predict a conjunction probability between the naval vessel and a further body in a body of water by implementing the method according to claim 10 or 11, anda control assembly (30) which is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the naval vessel a conjunction avoidance manoeuvre to be performed or to command a control system of the naval vessel to perform the conjunction avoidance manoeuvre, orwherein the control assembly (30) is adapted to compute, based on the conjunction prediction, and communicate to a pilot of the naval vessel a rendez-vous manoeuvre or to command a control system of the ground vehicle to perform the rendez-vous manoeuvre.

17. Method (1000-5000; 7000) of predicting a conjunction between an intercepting means (50H, 53) and a moving target (T) comprising the steps of:acquiring (2001; 2001 A, B,E, F— 2001 A,B,E,F) data defining at least one last known state and last known covariance of the intercepting means and of the target;computing (2003-2007; 2003A, B,E, F - 2007A, B,E, F) a plurality of Gaussian Mixture Elements (GMEs) based on the acquired data;computing (3001-3003; 3001 A, B, E, F - 3003A, B, E, F) a Gaussian Mixture Model (GMM) of an initial state based on the GMEs, the initial state being a state of the intercepting means and of the target at an initial time instant defining a lower boundary of a conjunction time interval comprising a time of closest approach of the intercepting means and of the target;computing (3005-3015; 3005A,B,E,F - 3015A,B,E,F) Taylor expansions of dynamic equations of motions of the intercepting means and of the target for each GME in the GMM;computing (3017-3025; 3017 A, B, E,F - 3025A, B,E, F) a final state and a GMM of the final state of the intercepting means and of the target based on the Taylor expansions, the final state being a state of the intercepting means and of the target at a final time instant defining a higher boundary of the conjunction time interval;computing (3039; 3039A, B, E, F) a conjunction probability rate as a function of time based on the initial and final states and GMMs of the initial and final states, andcomputing (3045; 3045A,B,E,F) a conjunction probability by integrating the conjunction probability rate over the conjunction time interval, anddetermining (1015; 7015) a possible conjunction if the conjunction probability is higher than a predetermined threshold, characterized in thatwherein the step of computing (3017-3025; 3017 A, B, E, F -3025 A,B,E,F) a final state and a GMM of the final state of the body and of the target comprises:defining (3019; 3019A, B,E,F) a first matrix and a second matrix based on the Taylor expansion of state dynamic equations, the first matrix comprising coefficients of the Taylor expansion and the second matrix comprising exponents of the Taylor expansion,defining (3021; 3021 A, B, E,F) a covariance matrix based on the GMM of the initial state, andcomputing (3023; 3023A, B, E,F) the GMM of the final state based on a combination of the first matrix, the second matrix and the covariance matrix.

18. System (1G; 1H) of conjunction prediction and for at least one intercepting means (50H) or a propelling assembly (50G) to propel an intercepting means (53) comprising:a conjunction prediction assembly (10; 10H), comprising at least a computing element, which is adapted to predict a conjunction probability between the intercepting means and a moving target by implementing the method according to claim 17, anda control assembly (30) which is adapted to compute, based on the conjunction prediction, and communicate to the intercepting means a conjunction avoidance manoeuvre to be performed or to command to the propelling assembly a gun laying adjustment to avoid the conjunction, orwherein the mission control assembly (30) is adapted to compute, based on the conjunction prediction, and communicate to the intercepting means a conjunction reaching manoeuvre to be performed or to command to the propelling assembly a gun laying adjustment to reach the conjunction with the target.

19. Method (1000G-7000G) of predicting a conjunction between an intercepting means (50H, 53) and a static target (T) comprising the steps of:acquiring (2001 G) data defining at least one last known state and last known covariance of the intercepting means (50H, 53) and a position of the target (T);computing (2003G-2007G) a plurality of Gaussian Mixture Elements (GMEs) of the intercepting means based on the acquired data;computing (3001G-3003G) a Gaussian Mixture Model (GMM) of an initial state based on the GMEs, the initial state being a state of the intercepting means at an initial time instant defining a lower boundary of a conjunction time interval comprising a time of closest approach of the intercepting means to the target;computing (3005G-3015G) Taylor expansions of dynamic equations of motions of the intercepting means for each GME in the GMM;computing (3017-3025) a final state and a GMM of the final state of the intercepting means based on the Taylor expansions, the final state being a state of the intercepting means at a final time instant defining a higher boundary of the conjunction time interval;computing (3039G) a conjunction probability rate as a function of time based on the initial and final states and GMMs of the initial and final states, andcomputing (3045G) a Circular Error Probable (CEP) by integrating the conjunction probability rate over the conjunction time interval, anddetermining (1015G; 7015G) a possible conjunction if the CEP is higher than a predetermined threshold, characterized in thatwherein the step of computing (3017G-3025G) a final state and a GMM of the final state of the body comprises: defining (3019G) a first matrix and a second matrix based on the Taylor expansion of state dynamic equations, the first matrix comprising coefficients of the Taylor expansion and the second matrix comprising exponents of the Taylor expansion,defining (3021 G) a covariance matrix based on the GMM of the initial state, andcomputing (3023G) the GMM of the final state based on a combination of the first matrix, the second matrix and the covariance matrix.

20. System (1G; 1H) of conjunction prediction and for at least one intercepting means (50H) or a propellingassembly (50G) to propel an intercepting means (53) comprising:a conjunction prediction assembly (10; 10H), comprising at least a computing element, which is adapted to predict a CEP between the intercepting means and a moving target by implementing the method according to claim 19, and a control assembly (30) which is adapted to compute, based on the conjunction prediction, and communicate to the intercepting means a target avoidance manoeuvre to be performed or to command to the propelling assembly a gun laying adjustment to avoid the target, orwherein the mission control assembly (30) is adapted to compute, based on the conjunction prediction, and communicate to the intercepting means a target hitting manoeuvre to be performed or to command to the propelling assembly a gun laying adjustment to hit the target.

Images