Determination of the Probability of Collision of Near-Earth Space Objects Taking into Account Form and Orientation

A method for assessing the probability of collision of two observed space objects in near-earth space is proposed. The solution is based on the method of statistical modeling, which made it possible to take into account the shape and orientation of approaching objects. The developed method provides for the automatic stop of calculations when the required accuracy and reliability of the solution is achieved. An analysis of existing approaches to assessing the probability of collision is given. A calculation was made for the real approach of the vehicles CHINASAT-20 and SHIJIAN-17.


INTRODUCTION
Operators of modern spacecraft (SCs) in the course of their activities are forced to take into account the threat of collision with objects in near-Earth space (NES). In order to provide early warning of dangerous situations in different countries, the APSOS (Roscosmos), CRASS (ESOC), SOCRATES (Analytical Graphics Inc.), STK Advanced CAT (Analytical Graphics Inc.), etc., specialized systems have been created. Every year, space-flight control centers register up to 70 potentially dangerous encounters with individual active spacecraft [1].
In addition to identifying the fact of an impending rapprochement, in order to make a decision on the implementation of an evasion maneuver, it is necessary to assess the degree of threat. The most appropriate criterion is the probability of collisions. Several methods are known to determine the probability under the following assumptions.
2) Within the limits of the approach interval, the movement is rectilinear.
3) The components of the SO velocities are determined without errors. 4) Errors in determining the motion parameters are distributed according to the normal law.
In [2], the problem is reduced to finding the probability of the relative position vector at the point of closest approach to the sphere of a given radius. The probability is determined based on the calculation of the volume integral using the approximation by cubes of small size, or using asymptotic series.
In [3,4], an analytical formula for calculating the probability is derived, while the general assumptions are supplemented by the following one.
5) The size of the scattering ellipsoid of an approaching SO significantly exceeds the dimensions of a protected space object (PSO).
This assumption prevents the method from being applied to large-sized PSOs, such as the ISS. In [5], an alternative approach to solving the problem is proposed, replacing assumption 5 with the following assumption.
6) The linear dimensions of the spacecraft being protected are much larger than the dimensions of the approaching objects.
Here, potentially dangerous SOs are approximated by material points. The study demonstrates that, for the moment of the closest approach in the direction of the relative velocity vector, the uncertainty of the relative position vector is equal to 0. This fact makes it possible to reduce the solution of the three-dimensional problem to the two-dimensional case by choosing a new coordinate system, one of the axes of which is collinear to the relative velocity vector.
The assumptions of the aforementioned approaches can significantly affect the value of the collision probability, as will be shown below, taking into account the real shape and orientation of the SO changes the calculation result several times. Assumption 2) does not allow one to evaluate the encounters, the duration of which takes a long time. The rest of the assumptions also make a significant contribution to the calculation error. These restrictions can be removed by resorting to statistical modeling methods (Monte Carlo). The use of the Monte Carlo method for calculating the probability of collision in the NES is found in several works, for example, in [6][7][8]; however, the algorithm for taking into account the shape and orientation of the SO is used only in study [9]. For this reason, the task was set as follows.

PROBLEM STATEMENT AND GENERAL APPROACH TO SOLUTION
Let the fact of convergence of two spacecraft be revealed on the time interval . At the moment of the beginning of the approach, state vectors X I (t 1 ) and X II (t 1 ), covariance matrices K I (t 1 ) and K II (t 1 ), the law in change of orientation, and the geometrical shapes of objects are known. It is necessary to find an estimate of collision probability satisfying a given level of accuracy and reliability (1) As was noted, it is possible to solve the problem without resorting to classical assumptions using the method of modeling random variables and statistical evaluation of their characteristics. In this case, the modeled random variable is the indicator of the event of a collision of two SOs: where is an elementary event and A is an object collision event.
The expected value event indicator is equal to the probability of the event. Therefore, in obtaining an estimate of the mathematical expectation ( ) based on a series of stochastic simulations of the approach process, an approximate value of the collision probability can be calculated: ≈ Required number of simulations n can be found based on inequality (1). Specific methods of calculation of n and estimates of the accuracy of the solution will be presented below.
The algorithm for calculating the probability of collision consists of the following steps: -generating an array of pseudorandom six-dimensional state vectors of the SO in accordance with and t K -revealing the fact of intersection of three-dimensional models of the SO so as to find potentially dangerous combinations on the interval and -determination of and evaluating the accuracy of the resulting value.
GENERATING AN ARRAY OF PSEUDORANDOM VECTORS At this stage, k realizations need to be obtained of the six-dimensional vector of the SO kinetic parameters distributed according to the normal law in accordance with mathematical expectation and covariance matrix K.
The first step is to generate a set of pseudorandom values uniformly distributed over the interval (U(0, 1)). The speed of the algorithm and the quality of the resulting random vectors depend on the choice of the uniform distribution generator. Taking into account the work on this topic [10,11,12], the Mersenne Twister (MT) algorithm was chosen. In Fig. 1 shows a uniform distribution U(0, 1) points in twodimensional space, obtained as a result of the operation of the software module for generating pseudorandom values based on the MT algorithm. Next, using the Box-Muller transform, a set of independent realizations of state vectors is obtained, distributed according to the normal law (0, 1) with zero mathematical expectation and unit covariance matrix. Figure 2 shows the result of the transition from U(0, 1) to (0, 1). The diagrams located at the top and on the right show the total number of points located in the corresponding range of coordinates. Based on these diagrams, you can get a visual representation of the distribution law.
To obtain vector taking into account the variance and interdependence of its components, it is necessary to carry out the following linear transformation: where is a diagonal matrix composed of the eigenvalues of covariance matrix K and is the matrix composed of the eigenvectors of covariance matrix K. Figure 3 shows the distribution of points of a possible position of a real spacecraft (CHINASAT-20) in projection onto the plane orbital coordinate system (CS). The final stage of generating random vectors is a parallel transfer according to a given mathematical expectation.

REVEALING THE INTERSECTION OF THREE-DIMENSIONAL MODELS OF SOs
Before detecting the fact of collision of objects of complex shapes for each pair of trajectories , it is advisable to select trajectories such that, when they are moved along, there is a con-  vergence at a distance less than where R I and R II are the radii of the spheres describing the real shape of the SO. Due to the fact that number of simulations within the assigned task can reach values of the order of 10 6 , the speed of this section of the algorithm is critical. Taking into account this fact, as well as the fact that objects in the approach interval are at a relatively small distance from each other and experience the same influence of disturbing factors, simplified mathematical models are used within the framework of preliminary filtering of pairs: -the linearized motion model for near-circular orbits [13]; and -the Keplerian model of motion for orbits with significant ellipticity.
As a result of elimination, there are many pairs of trajectories during the movement along which there is a collision of spherical objects. For this set, it is necessary to identify the fact of intersection of three-dimensional models on the interval The problem of creating a fast and efficient algorithm for detecting the intersection of two threedimensional shapes of a given shape is one of the central problems of robotics and computer games. As a result, at the moment, there are a number of proven methods and algorithms for solving this problem. The general approach is to split the calculation process into two phases: wide and narrow.

= +
In the wide phase, all objects are represented as bounding parallelepipeds, the sides of which are parallel to the coordinate axes. This simplification allows one to very quickly select objects of potential collision. The Sweep-and-Prune algorithm was chosen as the wide phase algorithm [14].
The Gilbert-Johnson-Keerthi (GJK) algorithm was chosen as the narrow-phase algorithm, which is widely used due to its stability, speed, and memory footprint [15]. The GJK algorithm uses mathematical concepts such as the Minkowski sum, support function, and simplex [16]. The operation of this algorithm can be briefly described as follows. A check is performed to see if the center of the CS belongs to an arbitrary initial simplex. If the check result is positive, then the figures intersect; otherwise, using the support function, an extreme point is calculated in the direction opposite to the radius vector of the point of the simplex closest to the center of the CS. Thereafter, the found point is added to the simplex and the point farthest from the center of the CS is removed. If the extreme point is one of the previously deleted or one of those included in the current simplex, then the computation ends.

ESTIMATION OF ACCURACY AND CONDITION OF AUTOMATIC
CALCULATION CESSATION The estimate of the mathematical expectation tends in probability to the true value when the number of simulations tends n to infinity. In this case, absolute  error the estimate obtained is proportional to As I A takes one of two possible values this random variable can be attributed to the class of Bernoulli numbers, for which, in turn, it is possible to directly calculate required number of realizations N* guaranteeing fulfillment of the inequality (1). There are several approaches to determining N* [17].
The roughest upper bound can be obtained based on the Chebyshev inequality, where X is a random value, is an expected value, and is variance. Let us select wherein Denoting we obtain that the inequality (2) is guaranteed to be fulfilled when For Bernoulli numbers, the variance does not exceed 0.25; therefore, the upper bound for the required number of simulations is (3) An alternative approach to determining N* is based on the central limit theorem: if independent random variables distributed according to the same law, then Based on this, we can write the approximate inequality ≈ which holds for where is the quantile of the normal distribution for the probability Using the upper bound for the variance of the Bernoulli numbers, we find that the inequality (2) is approximately fulfilled for (4) Comparing (4) and (3), we can conclude that, due to the absence of a confidence factor in denominator (4), magnitude grows much more slowly than On the other hand, it should not be forgotten that, in this case, inequality (2) is asymptotically true.
The third way of finding N* is based on Hefding's inequality for random variables lying in the range and, hence, N * can be found from the equation    The minimum value of N* is obtained by means of a method based on the conclusions of the central limit theorem. Even when using formula (4), the order of the number of simulations is 10 9 , which casts doubt on the applicability of the Monte Carlo method for calculating the collision probability. At the same time, it should be noted that, in formulas (3)-(5), it is assumed that the maximum possible value of the variance is In real calculations, the variance takes on a much smaller value and, as a consequence, the required accuracy is achieved with fewer simulations. This conclusion is confirmed by the graph in Fig. 4 showing the dependence of probability estimate on number of simulations N for several runs of the algorithm. The launches were carried out on the same initial data; a sphere was used as a three-dimensional model. Based on the graph, you can make an approximate conclusion about the rate of convergence of the method, an accuracy of is achieved at N = 1.6 × 10 6 .
On this basis, it is proposed to use the condition of automatic stopping of the calculation when the required accuracy and reliability of the solution is achieved. Starting with number of simulations N 0 , the current value of the variance is estimated and condition (4) is checked; if the condition is met, then the calculations are stopped.    Table 3 shows the results of calculating the probability of collision, taking into account the shape and orientation of the spacecraft. As part of testing, a generalized three-dimensional model of geostationary ±ε,     satellites was used, shown in Fig. 5. This model can be inscribed in a sphere with a radius of 15 m. Since solar panels, as a rule, have one degree of freedom, their shape was approximated by cylinders. It was assumed that the objects have a typical orientation for geostationary satellites: the axis of the solar panels is perpendicular to the orbital plane, and the main axis is directed to the center of the Earth. Similarly to the calculations for the spheres, three series of launches were performed for different values of required accuracy The scatter of the obtained estimates does not exceed the specified error, while the probability of collision of a complex-shaped SO was approximately five times less than the probability of collision of the describing spheres. Calculation time in single-threaded computation mode for ε = 2 × 10 -6 did not exceed 65 min.

CONCLUSIONS
The probability of collision of an SO is the main quantity influencing the decision to use an evasive maneuver. The safety of spaceflights and the rate of consumption of the spacecraft's working medium depend on the accuracy of the estimate of this value. Within the framework of the study, a method for calculating the probability was developed, which makes it possible to take into account the shape and orientation of approaching objects, as well as remove a number of classical assumptions. It was demonstrated that ignoring the real shape of objects significantly affects the accuracy of the estimate; for the test case, a fivefold discrepancy in the results was obtained. At the same time, when analyzing the convergence of spherical objects, the estimates of the developed method agree with the estimates of the classical method of Khutorovsky.

ε.
The proposed method provides an algorithm for automatically stopping calculations when a given accuracy and reliability of the solution is achieved. It should be noted that the required calculation time increases nonlinearly with a decrease in the value of the permissible error. In order to speed up the operation of the algorithm in the future, it is planned to replace the pseudorandom number generator of the Mersenne Twister with the Sobol sequence [18], which provides a higher rate of convergence of statistical modeling methods. Acceleration can also be achieved by resorting to a two-stage risk assessment, in which, at the first stage, the probability of collision of the describing spheres is calculated and, if the probability is higher than the threshold value, the actual shape and orientation of the approaching objects are taken into account. FUNDING This work was supported by the RUDN University Strategic Academic Leadership Program.