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.
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
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:
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.
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.
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.
| ||||||||||||||||||||||||||||||||||
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.
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.
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.
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)
(a)
(b)
(c)
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.
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.
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)
(a)
(b)
(c)
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.
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.
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.
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 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.
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).