International Journal of Aerospace Engineering
Volume 2019 (2019), Article ID 5316396, 23 pages
https://doi.org/10.1155/2019/5316396

Research Article

An Analytical Method for Calculating Short-Term Space Interception Probability

Wenzhe Ding,¹ Xinhong Li,² Hong Yang,¹ Zhibin Zhang,¹ and Jiping An¹

¹Graduate School, Space Engineering University, Beijing 101416, China
²Department of Aerospace Science and Technology, Space Engineering University, Beijing 101416, China

Correspondence should be addressed to Wenzhe Ding; dingwenzhe123@gmail.com

Received 9 July 2019; Revised 24 August 2019; Accepted 5 September 2019; Published 3 November 2019

Academic Editor: Paolo Gasbarri

Copyright © 2019 Wenzhe Ding et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

For small unguided space interceptors, the interception probability is an important index to evaluate their strike capability. However, the Monte Carlo (MC) based simulation method not only consumes a lot of computation time but also theoretically lacks interpretability. In this regard, this paper proposes an analytical method that can quickly and accurately calculate the short-term space interception probability. Firstly, by considering the effect of perturbation force on spacecraft as the effect of external force acceleration, the analytical calculation formula of the short-term state error covariance propagation of space target and interceptor is deduced. Next, by projecting the state error and rotating the coordinate system, the joint error distribution of the target and the interceptor in the calculation coordinate system at the time of closest approach (TCA) is obtained. Thereby convert the calculation of the space interception probability into the integral of the 2-dimensional probability density function in the circular domain. Then, the Laplace transform and Taylor expansion are used to obtain the exact power series expression and the maximum truncation error of the integral calculation, and the analytical calculation of the short-term space interception probability is realized. Finally, the effectiveness of the proposed method is verified by a simulation example. The proposed method can directly calculate the space interception probability according to the initial state error distribution of the target and interceptor, and the whole calculation process does not contain double integral operation. The proposed method has high computational efficiency, is suitable for on-orbit calculation, and provides effective support for the rapid evaluation of the strike capability of the space interceptor.

1. Introduction

The interceptor can be mounted on a satellite platform as a payload and to perform orbital interceptor missions such as space debris removal. Small nonguided interceptors (such as projectiles) have the advantages of lightweight and can be fired continuously, but due to the lack of guidance, it is necessary to make a reasonable assessment of their strike capability. As an important index to evaluate the striking ability of small unguided interceptors, interception probability plays an important role in the design, optimization, and performance evaluation of interceptors. Literature [1] systematically compares the calculation methods of circular error probability (CEP). Literature [2] applied the bootstrap method to resample small samples and calculated the cyclic error caused by repeated sampling, so as to analyze the influence of various error sources on the shooting accuracy of electromagnetic railgun (EMRG). In literature [3], the ground strike accuracy of EMRG is analyzed with Latin hypercube design (LHD). However, all the above researches are based on ground tests. Due to the unique motion of the interceptor in space, the above research methods are no longer applicable.

When calculating the probability of intercept in space, the initial error of interceptor should be analyzed first. The source of the interceptor’s initial error can be classified into two categories. One is the position error of interceptor caused by the position error of the satellite platform. The other is the interceptor velocity error caused by satellite velocity error, attitude adjustment error, launch velocity error, and other factors. Therefore, the initial error of the interceptor can be simplified into the state error of the interceptor after launching in orbit. Different from the ground strike test, the initial state of the target will have some initial errors due to the limited ground observation ability. Due to the influence of gravity and various kinds of perturbation, the initial states of the target and the interceptor will change continuously throughout the interception process. Therefore, it is necessary to study the distribution of state errors at TCA according to the initial state error distribution of the target and the interceptor. Then, when calculating the space interception probability, we first need to study the uncertainty propagation problem in orbit mechanics.

As for the uncertainty propagation problem in orbital mechanics, there is much relevant research literature [4–6]. For the linear propagation process of error, the literature [7] used linear covariance analysis to study the orbit error propagation. When the error propagation time is too long and the system has strong nonlinearity, the literature [8, 9] uses Monte Carlo to study the error propagation problem. In order to reduce the computational cost of MC, some related techniques for error propagation have been proposed, such as unscented transformation (UT) [10], polynomial chaos (PC) [11, 12], differential algebraic (DA) [13, 14], and Gaussian mixture model (GMM) [15–17]. For the short-term space interception task, the propagation time of the state errors of the target and the interceptor is short during the whole interception process. After the propagation, the orbital deviations of both are small, so the Tschauner-Hempel (TH) equation can be used to describe the relative motion between the “mean orbit” and the “deviation orbit” of the target and the interceptor, so as to obtain the analytical solution of the target and interceptor state error propagation.

After error propagation, target and interceptor reach the closest position at the TCA. At this point, there is a high relative speed between the target and the interceptor and the overall contact time between them is extremely short during the encounter. Therefore, short-term encounter model can be used to define the encounter process between the target and the interceptor [18]. The motion trajectories between the target and the interceptor in the encounter process are assumed to be straight lines, and the uncertainty of the velocities of the target and the interceptor in the encounter process is ignored. So, the probability calculation can be converted to solve the integral of the two-dimensional probability density function in the encounter plane. For this double integral, the literature [19] uses the Monte Carlo method for calculation. In order to speed up the calculation efficiency, Foster and Estes [20], Patera [21], and Alfano [22] proposed different numerical integration models, respectively. Besides, Chan [23] proposed an analytical method to approximate the integral formula. The method approximates the integral domain. When the probability density function (PDF) is close to the circular distribution, it can well approximate the actual probability value. However, when the ratio of the semiminor axis length of PDF elliptic distribution is large, the probability calculation will be biased. Compared with Chan’s method, the analytical calculation method based on the Laplace transform [24] avoids the approximation of the integral domain and is more suitable for the calculation of interception probability.

In summary, the calculation of space interception probability is a complex multistage problem. In this issue, uncertain propagation calculation of orbit and collision probability calculation is involved. Although an exact solution to the space interception probability can be obtained by a large number of Monte Carlo simulations, it often requires a large amount of computational cost and time cost. In this regard, this paper combines the characteristics of short-term space interception task to define “short-term” from two aspects: The first is that the overall time of the interception process is short. The second is that the overall contact time between the target and the interceptor is short. Based on the above definition, this paper presents an analytical method to calculate the short-time space interception probability. This method divides the space interception task into two stages: error propagation and instantaneous collision, and realizes the analytical calculation of the short-term space interception probability by modeling the two phases separately. In the error propagation phase, we establish the state error uncertainty propagation model for the target and interceptor based on the TH equation. By considering the perturbation effect as the external force acceleration, the analytical calculation formula of the uncertainty propagation is obtained. By formula, we can directly calculate the error covariance matrices of the state of the target and interceptor at the TCA. In the transient collision phase, we unify the state vectors and state errors of the target and interceptor at the TCA into the calculation coordinate system, thereby transforming the calculation of the interception probability from the 3-dimensional problem to the 2-dimensional problem. The Laplace transform and Taylor expansion are used to obtain the expression of the power series of the space interception probability, which realizes the accurate analytical calculation of the space interception probability. The main innovation of this paper is to provide a precise method to solve short-term space interception probability in a pure analytic way. The advantage of this method is that it combines the uncertainty propagation calculation of orbit with the space collision probability calculation. Through the multistage modeling, an accurate analytical method for calculating the short-term space interception probability is given. Using this method, the short-term space interception probability can be calculated directly according to the initial state distribution of the target and the interceptor. The Laplace transform method used in the modeling avoids the double integral operation required in the final probability calculation, which makes the method more suitable for on-orbit calculation and provides the effective support for the rapid evaluation of the interception capability of the space interceptor.

2. State Transition Matrix and Analytical Solution for Uncertainty Propagation Calculation

For the space interception probability problem, the error uncertainty propagation of target orbit and interceptor orbit should be considered first. In this section, the relative motion theory is used to study the error propagation of the two orbits. The mean orbit and deviation orbit of the spacecraft are regarded as two close space objects of “real” and “virtual,” respectively, so the propagation rule of orbital error can be described by relative motion theory.

2.1. Relative Motion Model and State Transition Matrix under the Elliptical Orbit

In this section, Local Vertical Local Horizontal (LVLH) coordinate system of spacecraft mean orbit is used to establish the relative motion equation. The origin of the coordinate system lies at the center of mass of the “real object.” The -axis points the center of the mass of the “real object” along the center of the earth. The -axis points to the direction of movement of the “real object” and is perpendicular to the -axis in the mean orbital plane of the spacecraft. The -axis satisfies the right-hand rule. Figure 1 shows the LVLH coordinate system of the spacecraft mean orbit.

Figure 1: The LVLH coordinate system of the spacecraft mean orbit.

In Figure 1, represents the distance between the “real object” and the center of the earth.

In the LVLH coordinate system of the spacecraft mean orbit, the relative motion of the “virtual object” and “real object” is

where represents the position-velocity relationship of the “virtual object” with respect to the “real object” in the LVLH coordinate system at time (the orbital deviation of the spacecraft), represents the external force acceleration, and and are the coefficient matrices, specifically expressed as follows:

where is the gravitational constant; is the orbital angular velocity of the “real object,” is the semimajor axis of the “real object orbit,” and is the orbital eccentricity of the “real object.”

In order to express the relative motion relationship between the space objects more simply, Tschauner and Hempel converted the independent variable of their relative motion state into the true anomaly of the “real object” and thus obtained the TH equation:

where and . The expressions of and are as follows:

where and [25], where is the orbital angular momentum of the “real object.”

According to the literature [26], the transition matrix from the relative state to is

The state transition matrix for is , where the expression of is as follows:

where , , , and .

So,

where is the state transition matrix in regard to .

2.2. Uncertainty Propagation of Errors

In the process of error uncertainty propagation, the influence of perturbation force on spacecraft is considered as the external acceleration. Since the main perturbations such as solar pressure and atmospheric resistance are modeled according to the white Gaussian noise process [27], the external acceleration in equation (1) can be assumed to be a random Gaussian process with mean and intensity :

Assuming that initial states of and are independent of each other and follow the Gaussian distribution in equation (9), then according to equation (1), the orbital error uncertainty of the spacecraft will be propagated by and .

where represents the mean of initial relative state distribution and is the covariance matrix of the initial relative state.

Further, the formula (1) is presented as a form of the stochastic differential equation:

where has the following properties:

Therefore, by solving equation (10), the relative state at time can be expressed as follows:

Since satisfies the Gaussian distribution and is a random Gaussian process, at any moment is a Gaussian random variable:

where and are distribution mean and covariance matrix of relative state at time , respectively.

According to equation (12), and can be divided into combinatorial forms, namely, and , where and are calculated according to the distribution propagation of , and and are calculated according to the process propagation of .

First, analyze the calculations of and . Since LVLH is a moving coordinate system, the orbital velocity deviation of spacecraft in equation (7) is the product of relative derivative, namely,

In the actual process, the measured orbital velocity deviation of the spacecraft corresponds to the absolute derivative in the Local Vertical Local Horizontal Inertial (LVLHI) coordinate system. Here, the LVLHI coordinate system is an inertial coordinate system, which coincides with the instantaneous LVLH coordinate system, but remains unchanged in the inertial space. Therefore, the orbital deviation of the spacecraft in LVLHI coordinate system can be expressed as follows:

and satisfy

So,

Suppose that the orbital deviation of the spacecraft measured at the initial moment satisfies the Gaussian distribution , then and can be calculated according to the following equation:

According to equations (11) and (12), and are, respectively,

Since the calculation of is complicated and contains integral operations, it is difficult to solve by using equation (19) directly. For this purpose, this paper completes the analytical calculation of by transforming the computational domain (convert in the time domain to in the true anomaly domain).

Similar to equation (11), we can get

So, , and we can get

According to equations (19) and (22), we can get

Therefore, we can first convert into according to equation (23), which is converted into the calculation of . And then, is obtained by state transformation. The analytical calculation method for will be given in Section 2.3.

2.3. Analytical Calculation of Error Uncertainty Propagation

As can be seen from Section 2.2, the calculation of uncertainty propagation of spacecraft orbit can be converted into the calculation of , , , and . It is generally assumed that the deviation means of the spacecraft orbit at the initial moment are zeros, then according to equation (18), we can get

From equations (7) and (18), we can get

The calculation for is converted to solve for firstly. Since , equation (21) can be converted as follows:

The calculation of requires the introduction of an adjoint system. According to the literature [28], we can get that

where is a constant matrix.

Therefore, the specific expressions of and can be obtained as follows:

By inverting the constant matrix , we can get

Then, substituting into equation (26), we can get

Further, simplify the integral equation in equation (31). Substituting equation (23) into equation (31), we can get

Since matrix contains the calculation of , direct calculation of equation (32) requires multiple integral operations. To simplify the calculation, the equation (32) is converted to the eccentric anomaly domain for operation.

According to the literature [25], it can be obtained that

where represents the eccentric anomaly of the “real object orbit” and is the eccentric anomaly of the “real object orbit” at the initial moment.

Substituting equation (33) into equation (32) to simplify integral operation, we can get

where is the integrand matrix, whose specific expression is given in the appendix [29].

Then, combining equations (22), (31), (32), and (34), can be directly calculated as follows:

Finally, according to equations (19), (24), (25), and (35), the spacecraft orbit error distribution at any time can be obtained directly by analytical calculation:

3. The Calculation Method of Space Interception Probability

According to the analytical calculation method for the uncertainty propagation of spacecraft orbits given in Section 2, the error covariance matrices of the target state and interceptor state at the TCA can be obtained. This section will use the states and state errors of the target and interceptor at the TCA to give a calculation method of the space interception probability. First, according to the interception task, we give the basic assumptions and define the encounter coordinate system. Next, a method of calculating the space interception probability is given in the calculation coordinate system by state error projection and coordinate system rotation.

3.1. Basic Assumptions and the Encounter Coordinate System

The following four calculation parameters are involved in the calculation of the space interception probability:

(1)The TCA of the target and the interceptor(2)The Miss Distance (MD) between the target and the interceptor at the TCA(3)Position and velocity vectors of target and interceptor at the TCA(4)Error covariance matrices of the target and interceptor at the TCA

For parameters 1, 2, and 3, they can be obtained in advance through the orbit extrapolation of the target and the interceptor according to the space interception mission. The 4th parameter can be calculated by error propagation according to the initial state error distribution of the target and the interceptor (see Section 2).

Since the contact process between the target and the interceptor is very complicated, the state error covariance matrices of the target and the interceptor will change during the contact process. Therefore, in order to ensure the accuracy of the calculation and simplify the calculation reasonably, we will make some reasonable assumptions about the contact process of the target and the interceptor according to the characteristics of the interception task.

In the interception mission, the relative speed between the target and the interceptor is very large (generally in km magnitude according to the different tasks) and the joint radius is small (generally in m level determined by the volumes of the target and interceptor). Therefore, the overall contact time between the target and the interceptor is very small. At the same time, ignoring the influence of the shape of the target and the interceptor, the following assumptions can be made:

(1)The motion of the target and the interceptor in the contact process can be assumed to be the linear motion with the constant velocity(2)The uncertainty of target and interceptor speed during the contact process can be ignored(3)The error ellipsoids of the target and the interceptor remain unchanged during the contact process(4)The target and the interceptor are spheres with the known radius

When the relative distance between the target and the interceptor is less than its joint radius, the interception is considered as a success.

Since the interception process is linearized, the error ellipsoid can be projected onto the two-dimensional plane for calculation. The position error covariance matrix at the TCA will be used in the calculation. Divide as shown in Figure 2, then

Figure 2: Partition of .

At the TCA, the relative position vector between the target distribution center and the interceptor distribution center is . According to the assumptions 1, 2, and 3, it can be obtained that the and the relative velocity are perpendicular to each other at the TCA. The plane that is perpendicular to and passes through can be regarded as the encounter plane. Establishing an encounter coordinate system in this plane can eliminate the error in the velocity direction so that the space interception probability can be calculated.

The encounter coordinate system takes the target as the origin, the -axis points to the projection point of the interceptor in the encounter plane, the -axis points to the direction, and the -axis satisfies the right-hand rule (as shown in Figure 3). The unit vectors of the three axes can be expressed as follows:

Figure 3: The encounter coordinate system.

Then, the transition matrix from the Earth Centered Inertial (ECI) coordinate system to the encounter coordinate system is

Let the transition matrix from the LVLHI coordinate system to the ECI coordinate system be , then the position vectors of the target and the interceptor in the encounter coordinate system are

where and are the position vectors of the target and the interceptor in the encounter coordinate system, respectively, and and are the position vectors of the target and the interceptor in the LVLHI coordinate system, respectively.

In Figure 3, is the encounter plane between the target and the interceptor.

3.2. Interception Probability Calculation Model

Through the calculation of error propagation in Section 2, the position error covariance matrices of target and interceptor in LVLHI coordinate system at the TCA can be obtained. After projecting the respective position errors of the target and the interceptor onto the encounter coordinate system, the position error covariance matrices can be expressed as follows:

where and are, respectively, the position error covariance matrices of the target and the interceptor in the LVLHI coordinate system at the TCA and and are the position error covariance matrices of target and interceptor in the encounter coordinate system at the TCA, respectively.

Since the position distribution of the target and the interceptor is independent of each other, the projected positions of the target and the interceptor in the encounter plane still satisfy the Gaussian distribution. So, the relative position vector between the target and the interceptor obeys the 2-dimensional Gaussian distribution in the encounter plane. Then, at the TCA, let the position vectors of the distribution center of the target and the interceptor in the encounter plane be and , respectively. So, the mean value of the relative position vector between the distribution centers of the target and interceptor in the encounter plane can be expressed as follows:

where is the relative distance between the distribution centers of the target and interceptor in the encounter plane.

Let the position error covariance matrices of the target and the interceptor in the encounter plane at the TCA be and , respectively. Then, refer to the composition of (as shown in Figure 2), the expressions of and are as follows:

where and represent the elements at the positions of the th row and the th column in the and , respectively, is equal to , and is equal to .

Then, at the TCA, the error covariance matrix of the relative position vector between the target and the interceptor in the encounter plane is

According to equations (42) and (44), the position error ellipsoids of the target and the interceptor at the TCA can be projected into the encounter plane to form a joint error elliptic domain. According to the assumption 4 in Section 3.1, set the target’s diameter and the interceptor’s diameter to and , respectively. The joint sphere composed of the target and the interceptor can be projected into the encounter plane to form a joint circular domain (the diameter of the circular domain is ). Figure 4 shows the joint error elliptic domain and the joint circular domain.

Figure 4: The joint error elliptic domain and the joint circular domain.

When the relative distance between the target and the interceptor in the encounter plane is less than the radius of the joint circle domain, the interception can be considered as a success. Therefore, it can be seen from Figure 4 that solving the space interception probability can be converted into solving the probability that the relative position vector of the target and the interceptor in the encounter plane is located in the joint circle domain.

According to equations (43) and (44), is a symmetric matrix. So, the joint error covariance matrix in the encounter plane at the TCA is

where is the correlation coefficient. are the standard deviations of the joint error ellipse along the -axis and -axis, respectively.

It can be seen from Figure 4 that when is not 0, the principal axis direction of the joint error ellipse domain is inconsistent with the orientation of the encounter coordinate system. Therefore, in order to facilitate the calculation, the encounter coordinate system needs to be rotated, so that the -axis of the rotating coordinate system is aligned with the long semiaxis of the error ellipse domain. The rotated coordinate system is called the calculation coordinate system . Figure 5 shows the rotation process.

Figure 5: Rotation of the encounter coordinate system to the calculation coordinate system.

In Figure 5, and are the expected values of the joint error ellipse along the -axis and -axis directions, respectively. is the rotation angle, and its calculation method is discussed in the following three cases:

(1)If and , then (2)If and , then (3)If , then is calculated by

After rotation, the joint error covariance matrix in the calculation coordinate system is

where and are the standard deviations of the joint error ellipse along with the -axis and -axis directions, respectively. At this time, the distribution parameters of the joint error in the calculation coordinate system are

Then, the 2-D PDF can be expressed as follows:

The interception probability can be expressed as the integral of the PDF in the joint circle domain:

4. Space Interception Probability Calculation Model Based on the Laplace Transform

In Section 3, the calculation of space interception probability has been transformed into solving the integral of 2-D PDF in the circle domain. For the calculation of the double integral, a suitable analytical method is needed to ensure the accuracy of the interception solution. In this regard, in view of the advantages of Laplace transform method for double integral calculation in reference [24], this section references literature [24] to give an analytical solution for the space intercept probability calculation.

4.1. Laplace Transform of the Integral Formula

First, the interception probability calculation formula represented by equation (50) is rewritten as a functional form of , where . The form of function is as follows:

Since each of the terms obtained by the function after the series expansion is similar in size and opposite in sign, it will cause the sum of the series term to become 0 for certain precision. At this point, the function should be set before the operation. Replacing with can avoid the elimination phenomenon. From the perspective of the Laplace transform, we use the form of function as , where is the regulatory factor. According to the literature [24, 30], the selectable factor is .

According to the Laplace transform theorem, the Laplace transform of function can be expressed as follows:

where is a complex variable.

According to the Fubini theorem [31], when , there is

And the final expression for can be

Define the function as follows:

Then, the Laplace transform of function can be expressed as follows:

where . Convert equation (56) to

where the expressions of , , , and are as follows:

Next, we perform the Taylor expansion on and achieve the Laplace inverse transformation on each term of the expansions. Through a series of calculations, the power series expression of the calculation formula of space interception probability can be obtained.

4.2. Power Series Expression of Space Interception Probability Calculation Formula

Let the value of be

Define , then the derivative of is

According to and equation (60), is differentiable in the complex plane. According to the definition of power series, can be expanded into the form of power series:

where

Since , the power series expression of is

Performing the Laplace inverse transformation on each term in equation (63), the power series form of the space interception probability calculation formula can be obtained:

Let , then (63) can be written as follows:

where

If the sum of the first terms of the power series is taken as the calculated value of the interception probability, there will be a truncation error. Let the sum of the first terms of the power series be , then the truncation error is

According to equation (67), when the number of terms of the power series and the diameter of the joint domain is given, is only related to . Define the variables and as follows:

According to equation (62), the following relationship can be obtained:

where the expression of is as follows:

Since , , , and are nonnegative, we can get

According to equation (71), we can get

Therefore, can be easily obtained by proof by contradiction and the value range of truncation error in the calculation of space interception probability is

If , then the following relations hold

Therefore, equation (73) can be rewritten as follows:

where is the maximum truncation error.

Therefore, by combining Sections 2 and 4 of this paper, we can obtain the entire calculation process of space interception probability, as shown in Algorithm 1.

1. Input:

    Orbital information and orbital error distribution of target and interceptor at the initial moment , ; Target and interceptor diameter , ; The number of items of power series N.

2. Error propagation calculation

    1) Get the TCA and MD from the orbit extrapolation according to the initial orbit information of the target and interceptor;

    2) Calculate the true anomaly of the target orbit and the true anomaly of the interceptor orbit corresponding to the TCA;

    3) Calculate and according to formula (25);

    4) Calculate the eccentric anomaly of the target orbit and the eccentric anomaly of the interceptor orbit corresponding to the TCA;

    5) Calculate , according to the appendix, formula (6), formula (30), and formula (35);

    6) According to equation (36), the state error covariance matrix and of the target and the interceptor at the TCA are obtained from and .

3. Error projection and coordinate system rotation

    1) According to equations (41), (43), and (44), calculate the joint error covariance matrix of the target and the interceptor at the TCA in the encounter plane;

    2) Calculate the rotation angle according to equation (46), and obtain the joint error distribution parameters , , , and in the calculation coordinate system according to equation (48).

4. Space interception probability calculation

    1) Calculate , , , and according to equation (58), and calculate ~ according to equation (66) and the term number N of power series;

    2) Calculate space interception probability and maximum truncation error according to equations (65) and (75).

5. Output:

    Space interception probability and maximum truncation error .

Algorithm 1: The complete process of analytical calculation of short-term space interception probability.

5. Simulation Example

In this section, a simulation example is given to illustrate the correctness of the proposed algorithm. Take the gravitational constant for , and the initial orbital parameters of the target and interceptor are given in Table 1.

Table 1: Initial orbital parameters of the target and interceptor.

Orbital parameter
Target orbit
Interceptor orbit

Semimajor axis
16545
36834.9

Eccentricity
0.4
0.676229

Orbit inclination
63
114.433

RAAN
25
191.606

Argument of perigee
20
111.057

Initial true anomaly
30
287.59

Select the integration step size to 0.1 s, and the interception process of the interceptor to the target can be obtained by the high precision orbit prediction (HPOP) model in the Satellite Tool Kit (STK). Figure 6 shows the entire interception process.

Figure 6: Interception process of interceptor to target in the J2000 coordinate system.

The TCA of the target and the interceptor is 1500 seconds, which can be obtained from orbit extrapolation. Tables 2 and 3, respectively, show the positions and velocities of the target and interceptor at this time.

Table 2: Location parameters of the target and interceptor at the TCA.

Location in the J2000 coordinate system
(km)
(km)
(km)

Target
-4291.771159
4390.710760
11366.515836

Interceptor
-4291.821845
4390.691209
11366.550420

Table 3: Velocity parameters of the target and interceptor at the TCA.

Velocity in the J2000 coordinate system
(km/s)
(km/s)
(km/s)

Target
-5.681314
-2.060952
1.042846

Interceptor
6.981309
1.460946
0.057147

At this point, the relative distance and velocity between the target and the interceptor are 64.4006 m and 13.18 km/s, respectively.

In order to verify the error propagation calculation method in Section 2, the comparative analysis is performed by Monte Carlo simulation. Suppose the initial error distribution of the target and the interceptor is

According to the literature [32], the intensity of is

Take a time step of 15 s, and Monte Carlo has a sample size of 1000. The HPOP orbit extrapolation model is used for shooting simulation. The position and velocity error distribution of the target and interceptor corresponding to each time step in the 1500 seconds in the three axes of the LVLHI coordinate system are calculated. The simulation results and the analytical calculation results of this paper are displayed in Figures 7 and 8.

(a)

(b)

(c)

Figure 7: Variation of the position error mean square error of the target and interceptor in the LVLHI coordinate system. (a) Position error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system. (b) Position error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system. (c) Position error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system.

(a)

(b)

(c)

Figure 8: Variation of the velocity error mean square error of the target and interceptor in the LVLHI coordinate system. (a) Velocity error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system. (b) Velocity error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system. (c) Velocity error mean square error propagation curves of the target and interceptor in the -axis direction in the LVLHI coordinate system.

It can be seen from Figures 7 and 8 that, at the 1500 seconds, the target state error mean square errors through the analytical calculation are 47.782 m, 35.552 m, 24.469 m, 29.127 mm/s, 14.899 mm/s, and 13.248 mm/s. The target state error mean square errors through the Monte Carlo simulation are 47.439 m, 35.053 m, 24.116 m, 28.886 mm/s, 14.852 mm/s, and 13.422 mm/s. The relative errors between the analytical calculation results and the Monte Carlo simulation results are, respectively, -0.72%, -1.42%, -1.46%, -0.83%, -0.32%, and 1.30%. Figure 9 shows the change of the relative error of the position error mean square error (REPEMSE) and the relative error of the velocity error mean square error (REVEMSE) of the target during the entire interception process.

Figure 9: The change of the REPEMSE and the REVEMSE of the target.

Similarly, it can be obtained that, at the 1500 seconds, the interceptor state error mean square errors through the analytical calculation are 163.041 m, 145.838 m, 143.022 m, 132.027 mm/s, 87.486 mm/s and 83.769 mm/s; The interceptor state error mean square errors through the Monte Carlo simulation are 165.590 m, 144.780 m, 141.210 m, 134.125 mm/s, 87.132 mm/s, and 83.002 mm/s; The relative errors between the analytical calculation results and the Monte Carlo simulation results are, respectively, 1.54%, -0.73%, -1.28%, 1.56%, -0.41%, and -0.92%. Figure 10 shows the change of the REPEMSE and the REVEMSE of the interceptor during the entire interception process.

Figure 10: The change of the REPEMSE and the REVEMSE of the interceptor.

It can be seen from Figures 9 and 10 that the theoretical calculation errors of the state error of the target and the interceptor do not exceed 2% in the whole interception process. The validity of the error propagation calculation method in Section 2 is verified. In order to show the calculation results more clearly, Figures 11 and 12 further show the confidence ellipses of the position error of the target and the interceptor at the 1500 seconds. The confidence ellipses contain 95% of the Monte Carlo calculation samples.

(a)

(b)

(c)

Figure 11: Contrast display of the confidence ellipse of the position error of the target at the TCA. (a) Confidence ellipse of the position error of the target on the XY plane at the 1500 seconds. (b) Confidence ellipse of the position error of the target on the XZ plane at the 1500 seconds. (c) Confidence ellipse of the position error of the target on the YZ plane at the 1500 seconds.

(a)

(b)

(c)

Figure 12: Contrast display of the confidence ellipse of the position error of the interceptor at the TCA. (a) Confidence ellipse of the position error of the interceptor on the XY plane at the 1500 seconds. (b) Confidence ellipse of the position error of the interceptor on the XZ plane at the 1500 seconds. (c) Confidence ellipse of the position error of the interceptor on the YZ plane at the 1500 seconds.

From Figures 11 and 12, it can be seen that, up to plot accuracy, the error confidence ellipse obtained by the analytical calculation at the TCA almost coincides with the error confidence ellipse obtained by the Monte Carlo simulation. Therefore, in the process of calculating the space interception probability, the error propagation analytical calculation method in Section 2 can be directly used to obtain the error covariance matrices of the target and the interceptor state at the TCA.

Next, the error projection and coordinate system rotation are performed on the obtained error covariance matrices of the target and the interceptor. First, according to equations (41), (43), and (44), the joint error covariance matrix

of the target and the interceptor in the encounter plane at the TCA can be calculated. Then, according to equation (46), the rotation angle of the coordinate system can be calculated. Finally, the joint error distribution parameters , , , and of the target and interceptor in the calculation coordinate system can be calculated by equation (48). Set the target diameter and interceptor diameter to 10 m and 0.2 m, respectively. Then, the joint error ellipse (95% confidence interval) and the joint circular domain of the target and the interceptor in the calculation coordinate system can be obtained. Figure 13 shows the joint error ellipse and joint circular domain in the calculation coordinate system.

Figure 13: The joint error ellipse and joint circular domain in the calculation coordinate system.

According to equation (58), parameters , , , and can be calculated. Take the power series terms number . According to equations (65) and (75), the space interception probability and maximum truncation error are calculated, respectively. Table 4 shows the calculation results.

Table 4: Space interception probability and maximum truncation error through the Laplace transform.

1

2

3

4

5

6

7

8

It can be seen from Table 4 that as the number of terms of the power series increases, the truncation error of the interception probability obtained by the Laplace calculation method decreases continuously. When is greater than 4, the space interception probability obtained by the analytical calculation remains substantially unchanged. Therefore, for the example in this section, the number of power series items can be taken as 4.

In order to verify the correctness of the Laplace calculation method, the results are compared with the Monte Carlo method and the Chan method [23]. Take the term number of the infinite series and the term number of the power series to be 4, and the simulation numbers of MC method to be 20 million. Table 5 gives the calculation results of the three methods. It needs to be explained that the Monte Carlo simulation used here is only for testing the double integral calculation method for solving the equation (50) but not for testing the whole space interception probability calculation method in this paper. During the simulation, 20 million sets of samples are randomly generated according to the joint error distribution of the target and the interceptor in the calculation coordinate system. Count the number of samples in the joint circle domain, so as to obtain the Monte Carlo simulation results of the interception probability calculation.

Table 5: Interception probability obtained by the three methods.

Monte Carlo

Chan

Laplace transform

It can be seen from Table 5 that the interception probabilities obtained by the three methods are very close. The computational error between the Laplace method and the Chan method is almost zero that verifies the validity of the Laplace algorithm. Next, the calculation results obtained by the Laplace method and the Chan method are further compared in Table 6 below.

Table 6: The relative error of Chan’s and Laplace’s methods with respect to MC and the calculation time spent by the Laplace method and Chan method.

Relative error with respect to MC method (%)
Computing time (ms)

Chan
0.02467
0.71333

Laplace transform
0.02442
0.67015

Table 6 shows the relative error of Chan’s and Laplace’s methods with respect to MC and the calculation time spent by the Laplace method and Chan method. It can be seen that the Laplace method is superior to the Chan method belonging to the analytical method in terms of accuracy and time.

Finally, in order to test the effectiveness of the proposed algorithm, a complete Monte Carlo simulation is performed. According to the state error distribution of the target and interceptor at the initial moment, 20 million sets of target and interceptor state samples are randomly generated. The orbit extrapolation is performed by the HPOP model to obtain the position coordinates of the target and the interceptor in the J2000 coordinate system at the 1500 seconds. Calculate the relative position vector between the target and the interceptor in each group, and count the number of groups in the 20 million group that satisfy the relative position vector magnitude is smaller than the joint radius between the target and the interceptor. Figure 14 shows the distribution of 20 million relative position vectors and the confidence error ellipsoid containing 99.73% of samples at the 1500 seconds.

Figure 14: Relative position vector distribution between target and interceptor in the J2000 coordinate system.

According to statistics, the number of relative position vectors smaller than the combined sphere radius in the 20 million groups is 10659. The interception probability is . The relative error between the analytical result and the Monte Carlo result is 1.011%, which proves the effectiveness of the proposed method.

6. Conclusion

In this paper, an analytical method for calculating the probability of space interception in short term is given by multistage modeling. Firstly, using the relative motion theory, the analytical formulas for the error propagation of the space target and the interceptor in short term are derived. Then, by unifying the state vectors and state errors of the target and the interceptor at the TCA into the calculation coordinate system, the calculation of space interception probability is transformed into solving the integral of 2-D PDF in the circle domain. Finally, the power series expression of space interception probability is obtained by the Laplace transform and Taylor expansion, which yields an accurate analytical calculation method for space interception probability.

Simulation examples show the following:

(1)In the short-term interception process, the orbital deviations of the target and the interceptor are relatively small. The accurate state distribution of the target and interceptor at the TCA can be obtained by using the analytical formula of uncertainty propagation derived from the TH equation(2)During the instantaneous collision, there is a high relative speed between the target and the interceptor. The short-term encounter model can accurately define the encountering process between the target and the interceptor(3)When the interception probability is calculated, the analytical calculation method based on the Laplace transform is used to transform the integral calculation into the complex plane, which avoids the approximation of the integral domain. Compared with the Chan method, the calculation accuracy is improved, which can meet the calculation requirements of the space intercept probability(4)The method proposed in this paper can calculate the short-term space interception probability only by using the initial state distribution and joint diameter of the target and the interceptor. And the calculation process does not contain double integral operation, leading to high computational efficiency

Appendix

A. The Expression of

Define the symmetric matrix as follows:

Then, the expression of the specific element in is as follows:

where

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported in part by the National Defense Science and Technology Innovation Special Zone Project (contract number 18-H863-01-ZT-002-055-01).

References

Y. Wang, G. Yang, D. Yan, Y. Wang, and X. Song, “Comprehensive assessment algorithm for calculating CEP of positioning accuracy,” Measurement, vol. 47, pp. 255–263, 2014.

J. Zhang and W. An, “Assessing circular error probable when the errors are elliptical normal,” Journal of Statistical Computation and Simulation, vol. 82, no. 4, pp. 565–586, 2012.

P. Ma, Y. Zhou, X. Shang, and M. Yang, “Firing accuracy evaluation of electromagnetic railgun based on multicriteria optimal Latin hypercube design,” IEEE Transactions on Plasma Science, vol. 45, no. 7, pp. 1503–1511, 2017.

Y.-z. Luo and Z. Yang, “A review of uncertainty propagation in orbital mechanics,” Progress in Aerospace Sciences, vol. 89, pp. 23–39, 2017.

Z. Yang, Y.-z. Luo, V. Lappas, and A. Tsourdos, “Nonlinear analytical uncertainty propagation for relative motion near J2-perturbed elliptic orbits,” Journal of Guidance, Control, and Dynamics, vol. 41, no. 4, pp. 888–903, 2018.

Z. Yang, Y.-Z. Luo, and J. Zhang, “Nonlinear semi-analytical uncertainty propagation of trajectory under impulsive maneuvers,” Astrodynamics, vol. 3, no. 1, pp. 61–77, 2019.

D. K. Geller, “Linear covariance techniques for orbital rendezvous analysis and autonomous onboard mission planning,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 6, pp. 1404–1414, 2006.

C. Sabol, K. Hill, K. Alfriend, and T. Sukut, “Nonlinear effects in the correlation of tracks and covariance propagation,” Acta Astronautica, vol. 84, pp. 69–80, 2013.

V. Vittaldev and R. P. Russell, “Space object collision probability via Monte Carlo on the graphics processing unit,” The Journal of the Astronautical Sciences, vol. 64, no. 3, pp. 285–309, 2017.

S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, vol. 92, no. 3, pp. 401–422, 2004.

B. A. Jones, A. Doostan, and G. H. Born, “Nonlinear propagation of orbit uncertainty using non-intrusive polynomial chaos,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 2, pp. 430–444, 2013.

B. A. Jones, N. Parrish, and A. Doostan, “Postmaneuver collision probability estimation using sparse polynomial chaos expansions,” Journal of Guidance, Control, and Dynamics, vol. 38, no. 8, pp. 1425–1437, 2015.

M. Valli, R. Armellin, P. di Lizia, and M. R. Lavagna, “Nonlinear mapping of uncertainties in celestial mechanics,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 1, pp. 48–63, 2013.

A. Wittig, P. di Lizia, R. Armellin, K. Makino, F. Bernelli-Zazzera, and M. Berz, “Propagation of large uncertainty sets in orbital dynamics by automatic domain splitting,” Celestial Mechanics and Dynamical Astronomy, vol. 122, no. 3, pp. 239–261, 2015.

K. J. DeMars, R. H. Bishop, and M. K. Jah, “Entropy-based approach for uncertainty propagation of nonlinear dynamical systems,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 4, pp. 1047–1057, 2013.

K. Vishwajeet, P. Singla, and M. Jah, “Nonlinear uncertainty propagation for perturbed two-body orbits,” Journal of Guidance, Control, and Dynamics, vol. 37, no. 5, pp. 1415–1425, 2014.

V. Vittaldev and R. P. Russell, “Space object collision probability using multidirectional gaussian mixture models,” Journal of Guidance, Control, and Dynamics, vol. 39, no. 9, pp. 2163–2169, 2016.

V. T. Coppola, “Including velocity uncertainty in the probability of collision between space objects,” in AAS/AIAA Spaceflight Mechanics Meeting, pp. 2012–2247, American Astronautical Soc, 2012.

S. Alfano, “Satellite conjunction Monte–Carlo analysis,” in Proceedings of AAS/AIAA Spaceflight Mechanics Meeting, pp. 2009–2233, American Astronautical Soc, 2009.

J. L. Foster and H. S. Estes, A Parametric Analysis of Orbital Debris Collision Probability and Maneuver Rate for Space Debris, no. article 25898, NASA JSC, 1992.

R. P. Patera, “General method for calculating satellite collision probability,” Journal of Guidance, Control, and Dynamics, vol. 24, no. 4, pp. 716–722, 2001.

S. Alfano, “A numerical implementation of spherical object collision probability,” Journal of Astronautical Sciences, vol. 53, no. 1, pp. 103–109, 2005.

F. K. Chan, Spacecraft Collision Probability, Aerospace Press, El Segundo, CA, 2008.

R. Serra, D. Arzelier, M. Joldes, J.-B. Lasserre, A. Rondepierre, and B. Salvy, “Fast and accurate computation of orbital collision probability for short-term encounters,” Journal of Guidance, Control, and Dynamics, vol. 39, no. 5, pp. 1009–1021, 2016.

D. A. Vallado, Fundamentals of Astrodynamics and Applications, vol. 12, Springer Science & Business Media, 2001.

K. Yamanaka and F. Ankersen, “New state transition matrix for relative motion on an arbitrary elliptical orbit,” Journal of Guidance, Control, and Dynamics, vol. 25, no. 1, pp. 60–66, 2002.

M. E. Campbell, “Collision monitoring within satellite clusters,” IEEE Transactions on Control Systems Technology, vol. 13, no. 1, pp. 42–55, 2005.

H.-C. Cho, S. Y. Park, and K. H. Choi, “Application of analytic solution in relative motion to spacecraft formation flying in elliptic orbit,” Journal of Astronomy and Space Sciences, vol. 25, no. 3, pp. 255–266, 2008.

S. Lee, H. Lyu, and I. Hwang, “Analytical uncertainty propagation for satellite relative motion along elliptic orbits,” Journal of Guidance, Control, and Dynamics, vol. 39, no. 7, pp. 1593–1601, 2016.

W. Gawronski, J. Müller, and M. Reinhard, “Reduced cancellation in the evaluation of entire functions and applications to the error function,” SIAM Journal on Numerical Analysis, vol. 45, no. 6, pp. 2564–2576, 2007.

G. Grosche, V. Ziegler, E. Zeidler, and D. Ziegler, Eds., Teubner-Taschenbuch der Mathematik, Springer-Verlag, 2013.

G. L. Slater, S. M. Byram, and T. W. Williams, “Collision avoidance for satellites in formation flight,” Journal of Guidance, Control, and Dynamics, vol. 29, no. 5, pp. 1140–1146, 2006.