Line Sampling Procedure for Extensive Planetary Protection Analysis

Exploration missions to other planets have to satisfy planetary protection requirements to limit the probability of impacts between mission-related objects and celestial bodies, with the goal of reducing the risk of contaminating them with biological material coming from Earth. The verification of these requirements can become a lengthy and computationally expensive task when addressed with common methods such as Monte Carlo simulations, as they involve analysing the interplanetary trajectories and the uncertainties associated with them for time spans up to 100 years, and estimating small probabilities with strict confidence levels. This paper presents novel improvements of the line sampling method, already introduced for the verification of planetary protection requirements as a way to estimate the impact probabilities more efficiently and with greater accuracy than achieved with standard Monte Carlo. These newly developed techniques are presented, with the aim of making the analysis with Line Sampling more effective, and providing more information about the distribution of impacts in the initial uncertainty distribution: an algorithm to identify the time intervals where most close approaches are clustered, and an algorithm to improve the determination of the main sampling direction and increase the accuracy of the probability estimation.


Introduction
During interplanetary missions, launcher stages and inactive spacecraft are often left along orbits that may come back to the Earth or reach other celestial bodies, with the risk of impacting and contaminating them. For this reason, planetary protection policies require to verify that the probability of contaminating a celestial body due to non-scheduled impacts from mission-related objects is below a given threshold for any interplanetary mission. The maximum acceptable impact probability and the confidence level associated with it depend on the class of mission under examination [1], spanning from 10 −3 up to 50-100 years after launch for a generic Mars mission, to 10 −6 for sample return missions or missions to the Jovian moons, with a minimum confidence level of 99%. The precise propagation of the orbital state of these objects is therefore of great importance to preserve planets and moons that are promising for the development of extraterrestrial microbial life, such as Mars or Europa among many.
Given the chaotic behaviour of the orbital propagation problem, fully non-linear techniques, such as Monte Carlo (MC) simulations, are preferred to tackle this kind of problem thanks to their general and flexible way of approaching collision and impact probability estimation. This is the most common approach for planetary protection analysis, which requires high accuracy to verify that the design of an interplanetary mission is compliant to the international guidelines [2][3][4][5]. In particular, Colombo et al. [2] and Jehn et al. [3] exploit the expression of the confidence interval by Wilson [6] to estimate in advance the number of MC runs necessary to grant a desired confidence level (as requested by planetary protection policy).
However, more efficient sampling methods may increase the precision of the probability estimate, or reduce the amount of simulations, and thus the computational cost. Since the probability levels to be verified are generally low, the analysis can require a large number of standard MC simulations, resulting in high computational cost and time: the number of simulations required to estimate the probability, indeed, increases as the expected probability decreases and the confidence level to be guaranteed increases. This becomes relevant especially in the case of planetary protection requirements, as low levels of impact probability for each class of mission must be guaranteed with strict confidence levels.
For this reason, the line sampling (LS) method was selected to increase the precision of the probability estimate, or reduce the amount of simulations, and thus the computational cost with respect to standard MC. LS is a MC-based approach for the estimation of small probabilities whose main feature is the reduction of the multi-dimensional integration problem across the uncertainty domain to many onedimensional integrals which are evaluated analytically. These integrals are defined along lines, all parallel to a reference direction, that are used to sample the initial distribution. The reference direction is determined so that it points toward a region of interest of the domain (that is, a subset of initial conditions that will lead to an event under study, such as a failure or an impact with a celestial body): if this is properly chosen, the method can considerably increase the accuracy of the solution or reduce the number of required system simulations with respect to a standard MC.
LS was originally developed for the reliability analysis of failure in complex structural systems [7,8], and later applied to risk estimation for orbital conjunction analysis in combination with differential algebra (DA) [9]. In previous works, the authors adapted and applied LS to the estimation of impact probability of small objects with major celestial bodies, mainly for the verification of planetary protection requirements for exploration missions, but including also Near Earth Asteroids [10,11]. In this work, this application of the LS method was further developed to improve its effectiveness in the cases under study.
This manuscript presents, in the first part, the LS algorithm and the theory behind it, studying the formulation that is given in the literature. This theoretical formulation was expanded further in order to better characterise the behaviour of the method according to the cases under analysis. Later, a series of algorithms based on LS are introduced with the aim of creating a numerical procedure capable of analysing complex cases: • an algorithm used to improve a first guess for the reference sampling direction to get it closer to an ideal case, thus increasing the accuracy of the analysis; • an algorithm to identify time windows where impact regions could be found by analysing the close approaches that are recorded during a preliminary sampling; • an iterative procedure to explore all possible impact regions and sample them if possible, with the goal of obtaining a wider overview of the distribution of impact regions within the initial uncertainty distribution, and estimating impact probability with higher accuracy and efficiency.
Two test cases regarding planetary protection analysis of interplanetary missions are provided with their results, as examples to demonstrate the performance of the presented techniques singularly or in the context of the overall procedure: the launcher upper stage of the Solar Orbiter mission, and orbiter of the Mars Sample Return mission. The manuscript is organised as follows: Sects. 2 and 3 will describe the LS algorithm and the theory behind it, introducing new developments along with it; Sect. 4 will present the additional techniques that were devised to improve the original LS method and apply it to a multi-event analysis; Sect. 5 will show the results of the application of the presented methods to different test cases, comparing them with the performance of standard MC simulations; finally, Sect. 6 will summarise the main results and present some of the objectives for continuing this work in the future.

Line Sampling
The LS method is a MC-based approach with the objective of estimating small probabilities (probabilities of failure that are generally in the range of 10 −7 to 10 −6 [7,8]). It consists of four steps: (1) the mapping of random samples from the physical coordinate space into a normalised standard space; (2) the determination of the reference 1 3 direction ; (3) the probing of the region of interest along the lines following the reference direction; (4) the estimation of the event probability.
In this work, the event of interest is an impact between the object being propagated (more specifically, a spacecraft or a launcher stage) and a celestial body. The impact event is defined as the condition when the distance between the spacecraft and any celestial body is below the critical radius of said body at any time during the propagation. Thus, the impact region in the initial uncertainty domain is defined as the subset of initial conditions that will lead to an impact within the selected propagation window.
A summary of each step of LS is provided hereafter, along with an introduction to the choices made for the numerical implementation of the LS technique.

Mapping onto the Standard Normal Space
The first step of the LS procedure is the definition of this transformation: all the random vectors x ∈ ℝ d of physical coordinates (position and velocity) that are drawn from the nominal uncertainty distribution in the following phases need to be mapped onto a normalised coordinate space.
Each vector is mapped to a new parameter vector ∈ ℝ d whose components j , j = 1, … , d are all associated with a unit Gaussian distribution, each with a probability density function (pdf) j defined as: The components of the vector have the property of being independent and identically distributed (iid) random variables, meaning that each has the same probability distribution as the others and all are mutually independent [8]. As a consequence, the joint pdf of these random parameters is equal to the product of all the single pdfs: Thanks to this property, this transformation grants efficiency to the method, especially for problems with high dimensionality: it reduces the problem to a series of one-dimensional analytical evaluations, thus enabling a simplification of the computation of the probability later in the procedure.
The direct and the inverse transformations, from the physical domain to the normalised one and vice versa, preserve the joint Cumulative Distribution Function (CDF) between the two coordinate spaces, and are defined as: The Journal of the Astronautical Sciences (2022) 69:1537-1572 with Φ and F being the CDF of the unit Gaussian distribution and the input uncertainty distribution of the problem, respectively.
Following the definition of the pdf , the joint CDF Φ is du is the error function. State and manoeuvre uncertainty is usually represented in Gaussian form in planetary protection analysis. For this reason, the Rosenblatt transformation [12], which is defined specifically for multivariate Normal distributions, is applied in this work: for Gaussian-distributed uncertainty parameters (as in the cases under study), both the direct and the inverse transformations (respectively Eqs. (4) and (5)) become linear [8,12]. The transformation is reported in Appendix A.

Determination of the Reference Direction
The reference direction can be determined in different ways [8]. In this work, it is determined as the direction of a normalised "centre of mass" of the impact region. This region is approximated by applying the Metropolis-Hastings algorithm [13,14] to generate a Monte Carlo Markov Chain (MCMC) lying entirely in the impact subdomain starting from an initial condition within it. MCMC simulation is a method for generating samples conditional on a region satisfying a given condition, according to any given probability distribution described by the pdf p( ) . The algorithm to generate a sequence of N S samples from a given sample u drawn from the distribution p( ) is briefly explained in [15]: 1. generate the sample by randomly sampling a user-defined "proposal" pdf p * ( u ); 2. compute the ratio r = p( )∕p( u ); 3. set ̃ = with probability min(1, r) and ̃ = u with the probability 1 − min(1, r) , where ̃ is the candidate for the next element of the chain; 4. check whether the candidate ̃ lies in the region of interest I or not: if ̃ ∈ I , accept it as the next sample u+1 =̃ ; else, reject it and take the current sample as the next sample u+1 = u .
The starting condition of the MCMC can also be found in different ways (e.g. numerical minimisation of the distance to the planet, or prior knowledge). In this work, the initial point is determined via a preliminary sampling, as will be described in Sect. 4. It is important to notice that the specific method used to determine the sampling direction is not a core phase of the LS method, thus representing one of the degrees of freedom in the implementation of the whole procedure.
After the impact region has been populated with N S samples, the reference direction is computed in the standard normal space as where u , u = 1, … , N S are the points of the Markov chain made of N S samples converted from the physical space into the standard normal space. The simulations performed for the Markov chain require additional computational effort with respect to standard MC methods. Nevertheless, this option provides a good coverage of the impact region and a resulting better accuracy of the final probability estimate.

Line Sampling
After determining the reference sampling direction, N T initial conditions x k , k = 1, … , N T are randomly drawn from the nominal uncertainty distribution and then mapped to standard normal coordinates as k using the transformation in Eq. (4). For each sample in the standard normal space, a line starting from k , k = 1, … , N T and parallel to is defined according to the parameter c k , such that where ⟨ , k ⟩ is the scalar product between and k . In this way, the problem is reduced to a series of one-dimensional evaluations associated with each sample, with c k also normally distributed in the standard space.
The standard domain is then explored along each line with the aim of identifying the impact region of an uncertainty set. For the problem under study, the occurrence of the events is expressed by the value of the non-dimensional performance function where Y is a function that maps the value of the parameter c (corresponding to an initial condition 0 ) to a performance index, and R imp is the selected critical distance. In this work, Y is the minimum distance from the celestial body of interest (e.g. the Earth) in the given time window, while R imp represents the radius around the planet that is selected to define an impact (not necessarily equal to the planet radius). According to this definition, it follows that The Journal of the Astronautical Sciences (2022) 69:1537-1572 The limit condition g ( (c)) = 0 is given by a trajectory reaching a minimum distance from the planet equal to the impact radius R imp .
The performance function is evaluated iteratively along each line k = c k + k ⟂ to identify the values of c k corresponding to its intersections with the impact region, as displayed in Fig. 1. Due to the nature of the problem under analysis (that is, a single close approach event within a given time interval), the hypothesis is made that a maximum of two intersections between each line and the impact region are found, meaning that two values of c k exist for each standard normal random sample k , k = 1, … , N T where the performance function is equal to zero: g (c k 1 ) = 0 and g (c k 2 ) = 0 . These two values represent states leading to the limit condition, while the values of c k in between represent states leading to impacts: The main hypothesis is valid when the impact region extends across the uncertainty domain and can be approximated as a flat or slightly curved surface. On the contrary, if the impact region is limited in size, the line being sampled may not intersect the impact region and no root c k is found. Finally, if the impact region has a twisted shape or disconnected shape due to many close approaches, more than two intersections (or none) with each line may be found.

Estimation of the Impact Probability
Once the values (c k 1 , c k 2 ), k = 1, … , N T are known for all the sampling lines, the unit Gaussian CDF provides each random initial condition k , k = 1, … , N T with the conditional impact probability P k (I) . Thanks to the properties of the normalised coordinate space, each P k (I) can be computed as a one-dimensional integrals along the corresponding line: If no intersections between a sampling line and the impact region were found, the value of the probability integral is equal to 0, and the associated conditional impact probability has a null contribution to the total probability.
The total probability of the event P (I) (which is identified with a planetary collision in the approach presented in this work) and the associated variance ̂2 P (I) are then approximated as follows: The total probability is estimated as the average of the conditional probabilities computed along each line (including those not intersecting the impact region as having a null partial probability). The variance is computed as the variance of all the conditional probabilities.
The number N T of sampling lines should be chosen to guarantee a desired level of accuracy, according to the equations given above. Possible criteria to guide this choice are discussed in Sect. 4.2.

Theoretical Basis
The LS method was introduced in previous works [11,16], where a general explanation of the theory behind it was presented, together with the results of its application to different test cases, to show how the choice of this sampling method can improve the efficiency of the MC simulations for planetary protection analysis.
In this work, the method is further analysed to understand how the shape of the impact region and the choice of the sampling direction influence the effectiveness of (12) The Journal of the Astronautical Sciences (2022) 69:1537-1572 LS with respect to standard MC. In the case of the LS, the literature already gives a qualitative estimation of its efficiency compared with the standard MC in terms of convergence rate [17]. A summary of it is reported in Sect. 3.1 to introduce the notation that will be used in Sect. 3.2.

Theoretical Formulation of the LS Method
The aim of this section is to give an introduction of the theoretical formulation of the LS method. In particular, it provides the analytical demonstration that LS has higher accuracy than standard MC. This treatise follows the explanation and the notation presented by Zio in [17], but is necessary in order to introduce the work presented in the next section.
In Sect. 3.2, the theory presented here will be expanded, in order to characterise further the performance of LS with respect to standard MC. In particular, an approximated formula will be developed, showing how the accuracy of LS is related to the choice of the reference sampling direction and to the shape of the impact region being sampled.
Let ∈ ℝ d be a random variable whose distribution is described by the multidimensional probability density function (pdf) q X ( ) , and let F be the subdomain of the variables leading to an event of interest (which can be seen as the failure of a system or, in this case, an impact with a celestial body). A performance func- The probability of the event F can be expressed as a multidimensional integral in the form where = (x 1 , … , x d ) T ∈ ℝ d is the vector of the uncertain variables of the system, and E[⋅] is defined as the expected value operator. Equation (15) shows that I F ( ) is also a random variable [17] and that P(F) can be interpreted as the expected value of I F ( ).
The variance 2 of I F ( ) is then defined as Since I 2 F ( ) = I F ( ) ∀ ∈ ℝ d due to I F ( ) being defined as a binary function, it follows that For the standard MC, an estimator P (F) of the probability P(F) as expressed in Eq. (15) is obtained by dividing the number of times that I F ( k ) = 1, k = 1, … , N by the total number of samples drawn N: If standard MC is interpreted as a Point Sampling method, in comparison with the LS, I F ( ) becomes the random variable that is sampled in order to estimate the probability of event F. In the application of LS, a coordinate transformation from the physical space to the standard normal space T X ∶ → brings as advantages the normalisation of the physical variables through the covariance matrix, and the possibility to express the multidimensional pdf as a product of d unit Gaussian standard distributions j ( j ): With reference to Fig. 1, in the d-dimensional standard normal space, the subdomain F is the subspace for which the samples = ( 1 , … , d ) T ∈ ℝ d−1 satisfy a given property (e.g. an impact with a planet or a system failure). With the assumption that with F 1 ∈ ℝ d−1 , in this way the region F corresponds to the values of such that the performance function g ( ) satisfies the relation Considering this change of variables and the definition in Eq. (20), the integral in Eq. (15) can be rewritten as and manipulated as follows: is the subset of the random variables which lead to a given result (e.g. an impact).
In case of the standard MC, Φ F 1 ( −1 ) is a discrete random variable equal to I F ( ) ; as a consequence: This is clear from Fig. 1, where the sampling procedure is represented highlighting the boundary of the region corresponding to the event F. As a consequence: The consequence of these properties is visible when considering the definition of variance of an estimator for the two methods. An estimator P (F) of the probability P(F) as expressed in Eq. (22) can be computed as where k = 1, … , N T are iid samples in the standard normal coordinate space. Given the generic definition of variance for P(F) following Eq. (22) as the variance of the estimator P (F) is defined as meaning that the variance of the estimator directly depends on the variance of the random variable Φ F 1 ( −1 ) : , the previous relation can be extended as follows: is the definition of the variance of the standard MC as given in Eq. (17), one can conclude from Eq. (27) that the variance obtained by the LS method is always smaller than the one given by standard MC, or at least equal to it.
can be defined as a measure of the efficiency of the sampling method, with lower values of meaning a higher efficiency of the method in converging to the exact value of the probability. Equation (27) demonstrates that the c.o.v. of estimator in Eq. (25) as given by the LS method is always smaller than the one given by the standard MC, implying that the convergence rate of the LS is always faster than, or as fast as, that of the standard MC.

New Developments
As introduced in the previous section, here the theory behind the LS method is taken further, to obtain more insight about it and about which parameters affect the the accuracy the method provides. The analytical formulas showing that the accuracy of the LS is always equal or higher to the accuracy given by standard MC were manipulated as part of the research work on LS. From this work, a formula based on approximations was obtained, showing how the accuracy of LS depends from important parameters such as the shape of the impact region and the determination of the sampling direction. This formula will be demonstrated and discussed. The analytical development presented here follows the notation used in the previous section and taken from [17].
While in the case of the standard MC Eq. (27) is easy to treat, since , with the first term defined as an "average" value of c( −1 ) (represented as a dashed blue line in Fig. 2) such that , and the second term as a variation with respect to this average value.
The hypothesis is made that c( −1 ) represents a small variation with respect to the average value ĉ , as in the case of a quasi rectilinear border of the region F orthogonal to the sampling direction . Under this hypothesis, the integral in Eq. (29) can be rewritten as this means that a new estimation for the worst covariance given by the LS method (nominally, from Eq. (27), equal to the one given by the standard MC) was obtained, which takes into account the probability level through the term (ĉ) , and the shape of the region F and the direction of sampling through the term Δc( −1 ).
When the approximation of small c( −1 ) is valid, meaning that the region F has a regular shape and is distributed across the initial uncertainty and the sampling direction is chosen properly so that it points toward it, and the probability level is low, the term 2 (ĉ) ⋅ E −1 c 2 ( −1 ) is also small, and we can say that the variance given by the LS is below a value f P(F), Δc( −1 ) such that 2 (P(F)) LS ≤ f P(F), Δc( −1 ) ≤ 2 (P(F)) MC , thus increasing the convergence rate of LS with respect to standard MC. On the contrary, when the approximation does not hold (that is in cases with high probability levels, non-optimal sampling (30) direction, or irregurarly shaped impact regions), f P(F), Δc( −1 ) grows toward the covariance level of the MC.

Correction of the Sampling Direction
As pointed out both in the available literature and in the previous considerations, the sampling direction is one of the key parameters determining the accuracy and efficiency of the LS method, where the ideal case is represented by the sampling direction being (almost) orthogonal to the boundary of the impact region.
While it is generally impossible to determine an optimal sampling direction a priori, Zio et al. [17] suggest a strategy to identify it as the one minimising the variance 2 P (F) N T of the probability estimator P (F) N T . This strategy consists in an optimisation search having the variance 2 P (F) N T as the objective function to be minimised, requiring the iterative evaluation of hundreds or thousands of possible solutions guess and 2N T or 3N T system model evaluations to be carried out to calculate the objective function 2 P (F) N T for each proposed solution. Therefore, the computational effort associated to this technique could be prohibitive with a system model code requiring hours or even days to run a single simulation.
The method that was developed for this work, instead, goes through a pre-processing phase used to correct an initial guess for the sampling direction guess . This guess solution is obtained from the set of N S samples generated in the initial Markov Chain already described in Sect. 2.2, and is then used in a short LS analysis over a low number of samples to gain information about the general position and shape of the impact region; the solutions given by this preliminary sampling are used to approximate the impact region as a hyperplane according to a multilinear regression (assuming this approximation is valid), thus using the norm vector orthogonal to the hyperplane as new sampling direction.
The full algorithm represented in Fig. 3 follows the following steps: 1. The initial Markov Chain is performed and a guess for the sampling direction guess is found; 2. A short LS is performed using a few initial samples (also drawn from the Markov Chain itself), obtaining the corresponding c k values at the boundary of the impact region; 3. The reference frame is rotated to an orthonormal base aligned with guess , in order to later define the hyperplane during the regression more easily; 4. The c k values are used to approximate the impact boundary of the impact region as a hyperplane, according to a multilinear regression scheme as where ∈ ℝ d is the vector collecting the coefficients defining the d-dimensional hyperplane and orthogonal to it; 5. The reference frame is rotated back to the initial frame of the standard normal coordinate space and the vector is set as the new sampling direction; 6. The main LS procedure is performed following the "corrected" direction.
The efficiency of this scheme relies on two main hypotheses: first, as already stated, that the impact region or its boundary can be approximated as a hyperplane; second, that the "corrected" sampling direction is "more optimal" than the initial one due to it being orthogonal to the boundary of the impact region, thus closer to the ideal case for the LS.
The application of this technique and its comparison with the LS using the standard procedure are shown in Sect. 5.

Definition of the Orthonormal Base
The orthonormal base used in the multilinear regression is generated according to the Gram-Schmidt orthogonalisation process [18] starting from the guess sampling direction, in order to obtain a base that is aligned with guess .
This algorithm constructs an orthonormal basis starting from a set of d lin- Since an initial set is needed, it is set 1 = guess in order to construct it starting from guess only; the other vectors 1 , … , d of the initial set are generated as i = ( i , … , d , 1 , … , i−1 ), i = 1, … , d in order to avoid singularities.
The orthogonalisation algorithm proceeds as follows: The Journal of the Astronautical Sciences (2022) 69:1537-1572 The operation ⟨ , ⟩ denotes the inner or scalar product between two vectors , ∈ ℝ d , with ⟨ , ⟩ ⟨ , ⟩ being the projection of onto . In this way, = [ 1 … d ] is defined as an orthonormal base, with 1 parallel to guess , making it possible to express the regression hyperplane in a more convenient way.

Algorithm for Multilinear Regression
Given a set of n dependent values y i , i = 1, … , n and a set of n independent variables i = (x i1 , … , x id ) T ∈ ℝ d , i = 1, … , n , the linear regression model assumes that there exists a relationship f ( i ) = y i ∀i = 1, … , n which can be modelled linearly as with ∈ ℝ d being an error variable representing the random noise due to the approximation.
Defining the following quantities the coefficients are found by imposing the minimisation of the quadratic quantity which gives as solution The vector of coefficients also contains the components of the defining vector orthogonal to the hyperplane. Thus, in the correction procedure described in Sect. 4.1, the new sampling direction is set parallel to the vector defined as (b 1 , … , b d , −1) T in order to point toward the hyperplane approximating the impact region.

Multi-event Analysis
As will be presented in Sect. 5, LS can be more efficient (to reach the same accuracy) than the standard MC in estimating the impact probability for a single event, meaning for a close encounter with a specific body in a specific narrow time interval. However, using the same procedure to analyse multiple events (more close approaches with more bodies over an extended period) can lead to less accurate results (for the same number of samples) than the standard MC: as (34) the sampling direction is computed to be orthogonal or nearly orthogonal to one impact region, using the same to sample different regions in the initial dispersion will increase the variance of the probability estimates corresponding to those impact regions, since that sampling direction will intersect them at a larger angle when compared with any direction computed specifically for those cases.
In order to maintain a high level of accuracy throughout the analysis, a repetitive process is presented, as represented in Fig. 4. The proposed solution makes use of repeated LS procedures based on different sampling directions: these directions are identified according to close approach windows, that is the time intervals where close approaches with planets occur, which can be visualised in the example in Fig. 5: in the picture, each dot represents a fly-by (crossing of the sphere of influence, or SOI) plotted according to the epoch and the minimum is reported as a dot according to the epoch and the distance, and different close approach windows are reported as thin rectangles; the colour identifies the order of sampling, starting from the intervals with the lowest distance distance, while the thin rectangles represent the time intervals used to look for impact regions; the colour identifies the order that will be followed for the sampling, starting from the interval with the lowest distance. For each of these windows, a Markov Chain is started, so that a sampling direction can be determined in case an impact region is present.
The full procedure, as it is described in the scheme in Fig. 4, follows these steps: 1. A preliminary MC is performed, using a low number of initial samples, and every close approach occurring during the propagations is stored in memory: this initial phase allows to obtain a first estimate of how many and where the impact regions might be; 2. The data about the fly-bys is analysed and the close approach (CA) windows are identified and sorted according to the minimum distance (with priority given to the the window leading to the closest CA): this allows to focus on the impact regions corresponding to the most relevant impact events (if any); 3. A Markov Chain is started to search for an impact region inside the first close approach window: this allows to populate the impact region (if any) and determine a first guess for the sampling direction; 4. At the end of the Markov Chain (a) In case impacts were detected, A preliminary LS is performed using a first guess for the sampling direction, which is then corrected using the algorithm introduced in Sect. 4.1; this allows to obtain a new sampling direction that is ideally closer to the optimal case for LS, which in turns allows to improve its accuracy; (b) A new set of random initial conditions is drawn from the distribution and a complete LS is performed as described in the previous sections to identify the impact region and compute the corresponding impact probability; (c) In case no impacts were detected, the next CA window in the priority order is selected, and the procedure goes back to point 4 to start a new Markov chain; 5. After all CA windows have been analysed, the impact probability is given by the weighted sum of all the partial impact probabilities as computed in Eq. (13).
Despite being computationally more complex and less memory efficient (due to the larger number of parameters involved and to the larger output size) than a single LS analysis, this method can ideally offer a complete overview of the impact regions inside the initial uncertainty distribution, providing additional information with respect to a normal planetary protection analysis, thus being able to be used directly in mission design.
Many parameters can affect the efficiency of this method in terms of number of propagations (compared with the standard MC) and its reliability in identifying correctly the impact regions (if any). In particular, the number of samples of preliminary MC (determined by confidence level), the length of the Markov Chains, the number of sampling lines, the tolerances used in the iterative process to determine the intersections between each sampling line and the impact region were identified as key parameters in these regards. For each of them, a trade off is required: while higher values would increase the accuracy of the determination of the impact regions and, consequently, of the probability estimate, the impact on the computational load may make the LS process less efficient than the standard MC; on the contrary, lower values would reduce substantially the number of orbital propagations needed to identify and sample the CA windows, but would not ensure the correct recognition of the impact regions.

Algorithm for Merging Close Approach Intervals
Given n intervals I i = [t i 1 , t i 2 ] with t i 1 < t i 2 (corresponding, respectively, to the SOI entry and exit epochs of each close approach in the preliminary analysis), algorithm to merge the close approach intervals proceeds as follow: 1. All intervals are sorted from earliest to latest starting epoch t i 1 (SOI entry) and put into a main stack; 2. The interval on top of the main stack I 1 = [t 1 1 , t 1 2 ] is compared with the following one I 2 = [t 2 1 , t 2 2 ] ; (a) If t 1 1 ≤ t 2 2 AND t 2 1 ≤ t 1 2 , the two intervals I 1 and I 2 overlap, and thus are merged into a new I = [t 1 1 , t 2 2 ] which becomes the interval on top of the stack I 1 ; (b) Else, I 1 is moved to an output stack and the comparison starts again from the interval on top of the main stack; 3. The comparisons are carried out until the main stack is empty.
The maximum complexity of this algorithm is of order O(N log N) (e.g., as used for sorting the minimum and maximum range values), with N being the number of recorded CAs. However, since the procedure is applied only once, it does not increase the computational load of the whole analysis in a relevant way.

Application and Results
Two test cases were selected to show how the proposed procedure works and its performance in identifying impact regions shaped differently with different probability levels: 1. The launcher upper stage of Solar Orbiter: • the data refers to an old option for a launch in October 2018, later discarded during the mission design (initial data in Table 1 taken from [2]); • the analysis focuses on the trajectory of the launcher upper stage following the injection into the interplanetary transfer orbit (aiming to a fly-by with Venus) and the separation from the spacecraft; • the planetary protection requirement applied in this case is the same used to protect Mars for a generic mission, with a probability level of 10 −4 .

The Mars Sample Return mission:
• the analysis focuses on the return trajectory from Mars after performing an Earth-avoidance manoeuvre; • the uncertainty distribution is built assuming errors of 1 m on all position components and 5 m/s on all velocity components within a 3 confidence level • the planetary protection requirement applied in this case is the same used to protect Earth for all sample return missions, with a probability level of 10 −6 .
For both cases, the initial state uncertainty is expressed as a 6x6 covariance matrix over the state only, and all initial data (epochs, states and uncertainties) is reported, and defined with respect to the EME2000 inertial reference frame (Earth's Mean Equator and Equinox at 12:00 Terrestrial Time on 1 January 2000): the x-axis is aligned with the mean equinox, the z-axis is aligned with the Earth's celestial North Pole, and the y-axis is rotated by 90 • East about the celestial equator [19].
The propagations are carried out in normalised units (with reference length and time equal to 1 AU and 1 solar year, respectively) using the adaptive Dormand-Prince Runge-Kutta scheme of 8 th order (DOP853), with absolute and relative tolerances both set to 10 −12 .
For both cases, the stopping tolerance of the iterative process for the LS is 10 −3 with a maximum of 12 propagations per line.
The comparison between the standard MC and the proposed approach based on LS is performed by analysing the following parameters: • the total number of orbital propagations N prop ; • for the LS only, the total number of sampling lines N lines ; • the impact probability estimate P (I); • the sample standard deviation ̂ of P (I); The overall number of propagations N P is selected as a measure of the computational burden of the methods, while the standard deviation ̂ is instead used as indicator of the accuracy of the result, with lower values corresponding to lower variability [8,9].
The total computational time is also reported and discussed among the results, however is not included the evaluations of the efficiency of the method due to the variability intrinsic to the choice of the machine and its workload. Table 1 reports the initial state, associated covariance matrix, and reference epoch for this mission case, which are the same ones reported in [2].

Launcher Upper Stage of Solar Orbiter
A preliminary MC simulation was performed to obtain information about the CA windows for the multi-event analysis. The number of runs was determined using the probability, confidence level, and maximum number of samples reported in Table 2: that value was obtained by arbitrarily scaling the probability set by the planetary protection requirement by two orders of magnitude. The table also shows the number of close approaches that were recorded and the number of CA windows that were identified via the merging algorithm introduced in Sect. 4.2.1, while Fig. 6 shows their distribution in time. Figure 6 shows the distribution in time of the CA windows found with the preliminary MC analysis, indicating each encounter with a planet as a dot according to its epoch and minimum distance from the encountered planet. The color scale indicates the sampling priority that is given to each CA window, starting from the one with the closest CAs. In particular, Fig. 7 shows a close-up of the CA window that was assigned the highest priority, which, in this case, is also the only CA window where an impact region was found via the sampling with a Markov Chain. This was an expected result, since, in this launch option, the upper stage injects the spacecraft into a trajectory aiming for a direct gravity assist with Venus in the first year of the mission [20]. In the case of the launcher upper stage of Solar Orbiter, the LS was applied both without and with the correction of the sampling direction. The solution without the correction is shown in Figs. 8 and 9, while the one where the correction was applied is shown in Figs. 10 and 11 : in both cases, the figures show the Δv distribution as grey dots and the impact region as red (solution of the standard MC simulation) and green dots (boundary found via LS), and the convergence of the solutions given by standard MC and LS, in terms of estimated value of the impact probability and its  Table 3 reports the numerical solutions in both cases, comparing them with the solution given by the standard MC simulation. The results reported in Table 3 show the benefit given by the application of LS over standard MC, as the proposed method is able to reach the same variance of the probability estimate given by the standard MC using a much lower number of propagations when the sampling direction is chosen properly: not only the number of sampling lines is reduced by almost 8 times, but the number of necessary propagations is reduced as well by almost 6 times. In fact, a choice of the sampling direction not orthogonal or nearly orthogonal to the impact region affects not only the efficiency of the LS (with a higher number of sampling lines necessary to obtain the same accuracy), but also the accuracy of the iterative process itself, as in average 9 propagations per line are needed against the 7 propagations per line in case the sampling direction is corrected to be "more optimal". In both cases, the number of lines is initially determined based on the value of the probability set by the planetary protection requirement, and it is then increased to reach the same accuracy given by the standard MC in the non corrected case, or reduced when that accuracy is reached in the corrected case.
The differences in the values of the estimated impact probability are due both to the determination of the sampling direction, as already pointed out, and to the accuracy of the iterative process used to identify the boundary of the impact region. This also confirms the theoretical considerations made in Sects. 3.1 and 3.2 about the importance of a proper determination of the sampling direction for a good accuracy of the LS solution.
It must be pointed out that the difference between the probability values estimated by the two methods is due to LS identifying only the main impact region, as seen in Figs. 8 and 10, and ignoring the isolated impacts found by standard MC. These may belong to extremely thin impact regions, which may be captured by LS with different values of simulation parameters (e.g. number of preliminary MC runs, length of Markov Chains), or may be outliers, which would require a different approach to be sampled. However, since the impact probability associated with the main impact region contributes the most to the overall estimation, the values computed via MC and LS do not differ substantially. Table 4, instead, shows the computational load given by the whole multi-event procedure compared with the one given by the standard MC, for the analysis with the correction of the sampling direction, highlighting the various phases of the LS-based procedure, which is split into: In this case, a maximum limit of 200 samples per Markov chain was set. It is clear that the Markov Chains are the most demanding phases with 66% of the total computational load of LS, since each CA window is sampled to search for an impact region, while, as expected, the preliminary MC phase does not increase significantly the computational time. Table 5 reports the initial state, associated covariance matrix, and reference epoch for this mission case: the uncertainty was arbitrarily set to 1 m in position and 5 m/s in velocity, uniformly over all the components. A preliminary MC simulation was performed to obtain information about the CA windows for the multi-event analysis. The number of runs was determined using the probability, confidence level, and maximum number of samples reported in Table 6.

3
The table also shows the number of close approaches that were recorded and the number of CA windows that were identified, while Fig. 12 shows their distribution in time. Figure 12 shows the distribution in time of the CA windows found with the preliminary MC analysis. In particular, Fig. 13 shows a close-up of the CA window that was assigned the highest priority, which, in this case, is also the only CA window where an impact region was found via the sampling with a Markov Chain,     given by standard MC and LS, in terms of estimated value of the impact probability and its associated standard deviation. Table 7 reports the numerical solutions in both cases, comparing them with the solution given by the standard MC simulation. Also in this case, the benefit of using LS is visible, as both the number of propagations and the variance of the probability are reduced by orders of magnitude with respect to the standard MC in identifying the main impact region. In particular, this test case is very close to the optimal case of the application of LS: not only the probability value to be estimated is very low, but also the impact region has a planar shape extended across the uncertainty distribution. This alone allows a higher efficiency even without a correction of the sampling direction, as the LS is able to reach the same accuracy of MC with a number of propagations orders of magnitude lower. In turn, the correction of the sampling direction allows a further reduction of the number of propagation required and a refinement of the estimation of the impact probability, obtaining a direction that is almost orthogonal to the boundary of the impact region (these observations follow the ones already done in [21]).
Also in this case, the computational efficiency of the whole multi-event analysis is compared with the standard MC simulation, through the data reported in Table 8, highlighting the various phases of the LS-based procedure with the correction of the sampling direction, as defined in the previous test case. In this case, a maximum limit of 1000 samples per Markov chain was set. Again, the Markov Chains sampling of each CA window is the most demanding phase, with almost 80% of the total computational load of LS. However, both in this case and in the previous one, the complete LS-based procedure reduces the computational time with respect to the MC analysis.
Finally, Table 9 tries to quantify the effect on the numerical efficiency of the angular offset of with respect to the corrected direction: the angular offset is compared to the ratios of propagations and of sampling lines in the corrected and non corrected cases, for both tests. It is possible to observe what has already been pointed out while commenting the two test cases in the previous paragraphs, that is the effect of both the correction of the sampling direction and of the shape of the impact region. In the case of the Solar Orbiter mission, the angular offset is relatively large, which is reflected in the difference in the number of propagations and sampling lines up to one order of magnitude. In the case of the Mars Sample Return mission, the differences are lower, due to both the more favourable sampling direction (much closer to the corrected version with respected to the other test case) and the better shape of the impact region (which allows a higher efficiency in general). A more thorough analysis of this aspect is left for a future work.

Conclusions
This manuscript focused on the improvement of the uncertainty sampling techniques currently used in planetary protection analysis. The LS method was presented as an alternative to standard MC simulations, with the goal of reducing the computational load of the statistical analysis by sampling the initial uncertainty distribution in a more efficient way. The LS method uses lines following a reference direction instead of points to sample the uncertainty space, identifying the boundaries of the impact regions inside of it and computing the impact probability using integrals along these lines where they intersect the impact regions. Since these integrals are evaluated analytically, the LS method can provide an estimation of the impact probability that has a lower variance (thus can be considered more accurate) than the one given by standard MC when the same number of samples are used, or reach a similar accuracy using fewer samples. While previous works adapted LS for single events, the method was improved by designing new algorithms to allow its application to cases where multiple impacts with planets are possible over long periods of time. The performance of LS was analysed in comparison with standard MC both numerically, by applying both methods to different mission cases, and analytically, by developing an approximated formula which highlights how the accuracy of LS depends on the geometry of the impact regions found within the initial uncertainty domain and the choice of the sampling direction.
In particular, two novel algorithms were devised to support the application of LS to cases where multiple close approaches with planets, and thus multiple impact regions, are possible, as commonly found in planetary protection analysis. One algorithm uses the information gained from a preliminary MC analysis to identify the time intervals where most close approaches are clustered and define close approach windows: this ensures that all potential impact regions, when present, are explored with LS. The second algorithm additionally uses an initial guess for the sampling direction and corrects it to make it closer to the ideal case: this greatly improves the accuracy of the method without increasing significantly the computational load, since it allows to sample the impact region in a more optimal way. The procedure shown here has been developed with the main objectives of making planetary protection analysis more accurate and efficient. The results presented in this work have proven that the LS method is more efficient and accurate of the standard MC Simulation method currently employed in planetary protection analysis, confirming the outcomes of previous studies on the subject. In all the test cases, the LS was able to reach a similar level of accuracy as the standard MC while using a lower number of propagations, or, viceversa, it reached a higher accuracy when performing the same amount of orbital propagations. Future work will focus on different aspects of the techniques introduced here.
In the first place, the theoretical formulation of LS will be expanded further, with the aim of obtaining an analytical expression of the variance of LS correlating the number of sampling lines and the estimated probability level. This will allow to estimate in advance the number of LS runs necessary to reach a given confidence level, as it is already possible for the standard MC by exploiting the expression of the confidence interval by Wilson [6].
Secondly, efforts will be made for improving the multi-event analysis, in particular the identification of the candidate impact regions and their exploration via Markov's Chain, with the goal of making the identification of regions of interest in the initial uncertainty more efficient. In particular, a general study aimed at quantifying the sensitivity of the proposed approach on the various parameters affecting its performance shall be done: the number of runs in the preliminary MC, the length of the Markov chains, the sampling direction, and the iterative method to sample the lines are some of the most important factors on which the LS techniques work.
On a last note, more advanced sampling techniques will be explored and applied to the estimation of impact probability for planetary protection analysis.

Appendix: Mapping Based on Multivariate Normal Distribution (Rosenblatt's Transformation)
This appendix presents the transformation used to map the physical coordinates to the standard normal coordinates and viceversa, as described in Sect. 2.1, starting from a multivariate Normal distribution. The explanation presented here follows the one given in [12], using a different notation.
The transformation described here is based on the mapping shown in Eq. (3), which preserves the CDF between the starting uncertainty and the standard normal distribution that characterises the normalised coordinate space used in the LS algorithm. The latter is represented by the pdf and CDF already expressed in Eqs. � with average ∈ ℝ k and covariance matrix ∈ ℝ d × ℝ d , the CDF (defined as F( ) = ℙ( ≤ ), where ∼ N( ; , ) ) cannot be obtained in closed form, but only estimated numerically, due to the multidimensionality of the problem.
For this reason, the transformation is treated making use of conditional probabilities, as In this way, the mapping in Eq. (3) conserving the CDF of the distribution can be rewritten in a more simple way: where Λ (r) is the determinant of the minor of Λ defined as [ ij], i, j = 1, … , r ≤ k , and Λ (r) ij is the cofactor of the element (ij) of the minor Λ (r) . This way, the direct and the inverse transformations become linear and assume the forms, respectively: