# An Analytical Method for Calculating Short-Term Space Interception Probability

Authors: Wenzhe Ding; Xinhong Li; Hong Yang; Zhibin Zhang; Jiping An
Journal: International Journal of Aerospace Engineering
 (2019)

---

## 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.


---

## Body

## 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.” Thex-axis points the center of the mass of the “real object” along the center of the earth. The y-axis points to the direction of movement of the “real object” and is perpendicular to the x-axis in the mean orbital plane of the spacecraft. The z-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 Figure1, R 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(1)Ẋt=AtXt+But,where Xt=Xt,Yt,Zt,Ẋt,Ẏt,ŻtT represents the position-velocity relationship of the “virtual object” with respect to the “real object” in the LVLH coordinate system at time t (the orbital deviation of the spacecraft), ut=uxt,uyt,uztT represents the external force acceleration, and At and B are the coefficient matrices, specifically expressed as follows:
(2)At=03×3I3×3A1tA2t,B=03×3I3×3,A1t=ω2t+2μR3ω̇t0−ω̇tω2t−μR3000−μR3,A2t=02ωt0−2ωt00000,where μ is the gravitational constant; ωt=aμ1−e2/R2 is the orbital angular velocity of the “real object,” a is the semimajor axis of the “real object orbit,” and e 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:(3)X~′f=A~fX~f+B~fuf,where X~′f=X~f,Y~f,Z~f,X~′f,Y~′f,Z~′fT and uf=uxf,uyf,uzfT. The expressions of At and B are as follows:
(4)A~f=03x3I3x3A~1fA~2,B~f=03x31ωf3/2I3x3,A~1f=3ρf0000000−1,A~2=020−200000,where ρf=1+ecosf and ωf=ρ2fμ2/h3 [25], where h is the orbital angular momentum of the “real object.”According to the literature [26], the transition matrix from the relative state X~f to Xt is
(5)Tf=h3/2μρI3×303×3μesinfh3/2I3×3μρh3/2I3×3.The state transition matrix forX~f is ψf,f0=ΦfΦ−1f0, where the expression of Φf is as follows:
(6)Φf=0−c0−s3esJ−201s1+1ρ0−c1+1ρ3ρ2J000cρ00sρ0−c′0−s′3es′J+sρ2002c−e02s31−2esJ000−sρ00cρ=H3×6H3×6′,where s=ρsinf, c=ρcosf, s′=cosf+ecos2f, and J=∫f0f1/ρ2τdτ.So,(7)Xt=Tfψf,f0T−1f0Xt0=Γt,t0Xt0,where Γt,t0 is the state transition matrix in regard to Xt.


### 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 ut in equation (1) can be assumed to be a random Gaussian process with mean 03×1 and intensity U3×3:
(8)Eut=03×1,EutuTt′=Uδt−t′.Assuming that initial states ofXt0 and ut 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 Xt0 and ut.
(9)Xt0∼ℕXt0;X0,C0,where X0 represents the mean of initial relative state distribution and C0 is the covariance matrix of the initial relative state.Further, the formula (1) is presented as a form of the stochastic differential equation:
(10)dXt=AtXtdt+Bdβt,where βt has the following properties:
(11)Eβt=03×1,Eβt2−βt1βt2−βt1T=∫t1t2Udt,t1<t2.Therefore, by solving equation (10), the relative state Xt at time t can be expressed as follows:
(12)Xt=Γt,t0Xt0+∫t0tΓt,τBdβτ.SinceX0 satisfies the Gaussian distribution and ut is a random Gaussian process, Xt at any moment is a Gaussian random variable:
(13)Xt∼ℕXt;mt,Ct,where mt and Ct are distribution mean and covariance matrix of relative state Xt at time t, respectively.According to equation (12), mt and Ct can be divided into combinatorial forms, namely, mt=m1t+m2t and Ct=C1t+C2t, where m1t and C1t are calculated according to the distribution propagation of X0, and m2t and C2t are calculated according to the process propagation of ut.First, analyze the calculations ofm1t and C1t. Since LVLH is a moving coordinate system, the orbital velocity deviation of spacecraft in equation (7) is the product of relative derivative, namely,
(14)Xt=Xt,Yt,Zt,Ẋt,Ẏt,ŻtT=r,δrδtT.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:(15)XLIt=Xt,Yt,Zt,ΔXt,ΔYt,ΔZtT=r,drdtT.dr/dt and δr/δt satisfy
(16)drdt=δrδt+ω×r=Ẋt−ωtYt,Ẏt+ωtXt,ŻtT.So,(17)XLIt=ML→LIωtXt=I3×303×3M21I3×3,M21=0−ωt0ωt00000.Suppose that the orbital deviation of the spacecraft measured at the initial moment satisfies the Gaussian distributionXLIt0∼ℕXLIt0;X0LI,C0LI, then m1t and C1t can be calculated according to the following equation:
(18)m1t=EML→LIωtΓt,t0MLI→Lωt0XLIt0=ML→LIωtΓt,t0MLI→Lωt0EXLIt0=φt,t0X0LI,C1t=Eφt,t0XLIt0−m1tφt,t0XLIt0−m1tT=φt,t0EXLIt0XLIt0TφTt,t0−m1tm1Tt=φt,t0C0LIφTt,t0−m1tm1Tt.According to equations (11) and (12), m2t and C2t are, respectively,
(19)m2t=06×1,C2t=∫t0tΓt,τBUBTΓTt,τdτ.Since the calculation ofΓt,τ is complicated and contains integral operations, it is difficult to solve C2t by using equation (19) directly. For this purpose, this paper completes the analytical calculation of C2t by transforming the computational domain (convert ut in the time domain to uf in the true anomaly domain).Similar to equation (11), we can get
(20)Eβ~f=03×1,Eβ~f2−β~f1β~f2−β~f1T=∫f1f2U~fdf,f1<f2,(21)C~2f=∫f0fψf,σB~σU~σB~TσψTf,σdσ.So,C2t=TfC~2fTTf, and we can get
(22)C2t=Tf∫f0fψf,σB~σU~σB~TσψTf,σdσTTf=TfT−1f∫f0fΓtf,τσTσB~σU~σ⋅B~TσTTσΓTtf,τσdσT−1fTTTf=∫f0fΓtf,τσTσB~σU~σB~TσTTσΓTtf,τσdσ=∫t0tΓt,τ1ωστBU~στ1ωστBTΓTt,τωστdτ=∫t0tΓt,τ1ωστBU~στBTΓTt,τdτ.According to equations (19) and (22), we can get
(23)U~f=ωfU.Therefore, we can first convertut into uf according to equation (23), which is converted into the calculation of C~2f. And then, C2t is obtained by state transformation. The analytical calculation method for C~2f will be given in Section 2.3.


### 2.3. Analytical Calculation of Error Uncertainty Propagation

As can be seen from Section2.2, the calculation of uncertainty propagation of spacecraft orbit can be converted into the calculation of m1t, C1t, m2t, and C2t. 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
(24)m1t=06×1,C1t=φt,t0C0LIφTt,t0.From equations (7) and (18), we can get
(25)C1t=ML→LIωtTfψf,f0T−1f0MLI→Lωt0C0LI⋅ML→LIωtTfψf,f0T−1f0MLI→Lωt0T.The calculation forC2t is converted to solve for C~2f firstly. Since ψf,f0=ΦfΦ−1f0, equation (21) can be converted as follows:
(26)C~2f=Φf∫f0fΦ−1σB~σU~σB~TσΦ−1σTdσΦTf.The calculation ofΦ−1f requires the introduction of an adjoint system. According to the literature [28], we can get that
(27)G=γTfΦf,where G is a constant matrix.
(28)γf=A~2H3×6−H3×6′H3×6.Therefore, the specific expressions ofγf and G can be obtained as follows:
(29)γf=23sρ0−3cρ−e3J3ρ−1+e2−3esρ200e001000sρ00−cρ0−c0−s3esJ−201s1+1ρ0−c1+1ρ3ρ2J000cρ00sρ,G=0−e00−10e001000000010−100−e0100e0000−1000.By inverting the constant matrixG, we can get
(30)G−1=11−e20−e0010e00−10000000e2−10100−e0−100e00001−e2000.Then, substitutingΦ−1f=G−1γTf into equation (26), we can get
(31)C~2f=ΦfG−1∫f0fγTσB~σU~σ⋅B~TσγσdσG−1TΦTf=ΦfG−1∫f0f1ω3σHTσU~σHσdσG−1TΦTf.Further, simplify the integral equation in equation (31). Substituting equation (23) into equation (31), we can get
(32)∫f0f1ω3σHTσU~σHσdσ=h6μ4∫f0f1ρ4σHTσUHσdσ.Since matrixHf contains the calculation of J, 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
(33)df=1−e21/21−ecosEdE,cosf=cosE−e1−ecosE,sinf=1−e21/2sinE1−ecosE,ρf=1−e21−ecosE,J=1−e2−3/2E−esinE−E0+esinE0,where E represents the eccentric anomaly of the “real object orbit” and E0 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
(34)∫f0f1ρ4σHTσUHσdσ=∫E0EKεdε,where KE is the integrand matrix, whose specific expression is given in the appendix [29].Then, combining equations (22), (31), (32), and (34), C2t can be directly calculated as follows:
(35)C2t=h6μ4TfΦfG−1∫E0EKεdεG−1TΦTfTTf.Finally, according to equations (19), (24), (25), and (35), the spacecraft orbit error distribution at any time can be obtained directly by analytical calculation:
(36)mt=m1t+m2t=06×1,Ct=C1t+C2t=φt,t0C0LIφTt,t0+h6μ4TfΦfG−1∫E0EKεdεG−1TΦTfTTf.


## 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.” Thex-axis points the center of the mass of the “real object” along the center of the earth. The y-axis points to the direction of movement of the “real object” and is perpendicular to the x-axis in the mean orbital plane of the spacecraft. The z-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 Figure1, R 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(1)Ẋt=AtXt+But,where Xt=Xt,Yt,Zt,Ẋt,Ẏt,ŻtT represents the position-velocity relationship of the “virtual object” with respect to the “real object” in the LVLH coordinate system at time t (the orbital deviation of the spacecraft), ut=uxt,uyt,uztT represents the external force acceleration, and At and B are the coefficient matrices, specifically expressed as follows:
(2)At=03×3I3×3A1tA2t,B=03×3I3×3,A1t=ω2t+2μR3ω̇t0−ω̇tω2t−μR3000−μR3,A2t=02ωt0−2ωt00000,where μ is the gravitational constant; ωt=aμ1−e2/R2 is the orbital angular velocity of the “real object,” a is the semimajor axis of the “real object orbit,” and e 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:(3)X~′f=A~fX~f+B~fuf,where X~′f=X~f,Y~f,Z~f,X~′f,Y~′f,Z~′fT and uf=uxf,uyf,uzfT. The expressions of At and B are as follows:
(4)A~f=03x3I3x3A~1fA~2,B~f=03x31ωf3/2I3x3,A~1f=3ρf0000000−1,A~2=020−200000,where ρf=1+ecosf and ωf=ρ2fμ2/h3 [25], where h is the orbital angular momentum of the “real object.”According to the literature [26], the transition matrix from the relative state X~f to Xt is
(5)Tf=h3/2μρI3×303×3μesinfh3/2I3×3μρh3/2I3×3.The state transition matrix forX~f is ψf,f0=ΦfΦ−1f0, where the expression of Φf is as follows:
(6)Φf=0−c0−s3esJ−201s1+1ρ0−c1+1ρ3ρ2J000cρ00sρ0−c′0−s′3es′J+sρ2002c−e02s31−2esJ000−sρ00cρ=H3×6H3×6′,where s=ρsinf, c=ρcosf, s′=cosf+ecos2f, and J=∫f0f1/ρ2τdτ.So,(7)Xt=Tfψf,f0T−1f0Xt0=Γt,t0Xt0,where Γt,t0 is the state transition matrix in regard to Xt.


## 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 ut in equation (1) can be assumed to be a random Gaussian process with mean 03×1 and intensity U3×3:
(8)Eut=03×1,EutuTt′=Uδt−t′.Assuming that initial states ofXt0 and ut 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 Xt0 and ut.
(9)Xt0∼ℕXt0;X0,C0,where X0 represents the mean of initial relative state distribution and C0 is the covariance matrix of the initial relative state.Further, the formula (1) is presented as a form of the stochastic differential equation:
(10)dXt=AtXtdt+Bdβt,where βt has the following properties:
(11)Eβt=03×1,Eβt2−βt1βt2−βt1T=∫t1t2Udt,t1<t2.Therefore, by solving equation (10), the relative state Xt at time t can be expressed as follows:
(12)Xt=Γt,t0Xt0+∫t0tΓt,τBdβτ.SinceX0 satisfies the Gaussian distribution and ut is a random Gaussian process, Xt at any moment is a Gaussian random variable:
(13)Xt∼ℕXt;mt,Ct,where mt and Ct are distribution mean and covariance matrix of relative state Xt at time t, respectively.According to equation (12), mt and Ct can be divided into combinatorial forms, namely, mt=m1t+m2t and Ct=C1t+C2t, where m1t and C1t are calculated according to the distribution propagation of X0, and m2t and C2t are calculated according to the process propagation of ut.First, analyze the calculations ofm1t and C1t. Since LVLH is a moving coordinate system, the orbital velocity deviation of spacecraft in equation (7) is the product of relative derivative, namely,
(14)Xt=Xt,Yt,Zt,Ẋt,Ẏt,ŻtT=r,δrδtT.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:(15)XLIt=Xt,Yt,Zt,ΔXt,ΔYt,ΔZtT=r,drdtT.dr/dt and δr/δt satisfy
(16)drdt=δrδt+ω×r=Ẋt−ωtYt,Ẏt+ωtXt,ŻtT.So,(17)XLIt=ML→LIωtXt=I3×303×3M21I3×3,M21=0−ωt0ωt00000.Suppose that the orbital deviation of the spacecraft measured at the initial moment satisfies the Gaussian distributionXLIt0∼ℕXLIt0;X0LI,C0LI, then m1t and C1t can be calculated according to the following equation:
(18)m1t=EML→LIωtΓt,t0MLI→Lωt0XLIt0=ML→LIωtΓt,t0MLI→Lωt0EXLIt0=φt,t0X0LI,C1t=Eφt,t0XLIt0−m1tφt,t0XLIt0−m1tT=φt,t0EXLIt0XLIt0TφTt,t0−m1tm1Tt=φt,t0C0LIφTt,t0−m1tm1Tt.According to equations (11) and (12), m2t and C2t are, respectively,
(19)m2t=06×1,C2t=∫t0tΓt,τBUBTΓTt,τdτ.Since the calculation ofΓt,τ is complicated and contains integral operations, it is difficult to solve C2t by using equation (19) directly. For this purpose, this paper completes the analytical calculation of C2t by transforming the computational domain (convert ut in the time domain to uf in the true anomaly domain).Similar to equation (11), we can get
(20)Eβ~f=03×1,Eβ~f2−β~f1β~f2−β~f1T=∫f1f2U~fdf,f1<f2,(21)C~2f=∫f0fψf,σB~σU~σB~TσψTf,σdσ.So,C2t=TfC~2fTTf, and we can get
(22)C2t=Tf∫f0fψf,σB~σU~σB~TσψTf,σdσTTf=TfT−1f∫f0fΓtf,τσTσB~σU~σ⋅B~TσTTσΓTtf,τσdσT−1fTTTf=∫f0fΓtf,τσTσB~σU~σB~TσTTσΓTtf,τσdσ=∫t0tΓt,τ1ωστBU~στ1ωστBTΓTt,τωστdτ=∫t0tΓt,τ1ωστBU~στBTΓTt,τdτ.According to equations (19) and (22), we can get
(23)U~f=ωfU.Therefore, we can first convertut into uf according to equation (23), which is converted into the calculation of C~2f. And then, C2t is obtained by state transformation. The analytical calculation method for C~2f will be given in Section 2.3.


## 2.3. Analytical Calculation of Error Uncertainty Propagation

As can be seen from Section2.2, the calculation of uncertainty propagation of spacecraft orbit can be converted into the calculation of m1t, C1t, m2t, and C2t. 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
(24)m1t=06×1,C1t=φt,t0C0LIφTt,t0.From equations (7) and (18), we can get
(25)C1t=ML→LIωtTfψf,f0T−1f0MLI→Lωt0C0LI⋅ML→LIωtTfψf,f0T−1f0MLI→Lωt0T.The calculation forC2t is converted to solve for C~2f firstly. Since ψf,f0=ΦfΦ−1f0, equation (21) can be converted as follows:
(26)C~2f=Φf∫f0fΦ−1σB~σU~σB~TσΦ−1σTdσΦTf.The calculation ofΦ−1f requires the introduction of an adjoint system. According to the literature [28], we can get that
(27)G=γTfΦf,where G is a constant matrix.
(28)γf=A~2H3×6−H3×6′H3×6.Therefore, the specific expressions ofγf and G can be obtained as follows:
(29)γf=23sρ0−3cρ−e3J3ρ−1+e2−3esρ200e001000sρ00−cρ0−c0−s3esJ−201s1+1ρ0−c1+1ρ3ρ2J000cρ00sρ,G=0−e00−10e001000000010−100−e0100e0000−1000.By inverting the constant matrixG, we can get
(30)G−1=11−e20−e0010e00−10000000e2−10100−e0−100e00001−e2000.Then, substitutingΦ−1f=G−1γTf into equation (26), we can get
(31)C~2f=ΦfG−1∫f0fγTσB~σU~σ⋅B~TσγσdσG−1TΦTf=ΦfG−1∫f0f1ω3σHTσU~σHσdσG−1TΦTf.Further, simplify the integral equation in equation (31). Substituting equation (23) into equation (31), we can get
(32)∫f0f1ω3σHTσU~σHσdσ=h6μ4∫f0f1ρ4σHTσUHσdσ.Since matrixHf contains the calculation of J, 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
(33)df=1−e21/21−ecosEdE,cosf=cosE−e1−ecosE,sinf=1−e21/2sinE1−ecosE,ρf=1−e21−ecosE,J=1−e2−3/2E−esinE−E0+esinE0,where E represents the eccentric anomaly of the “real object orbit” and E0 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
(34)∫f0f1ρ4σHTσUHσdσ=∫E0EKεdε,where KE is the integrand matrix, whose specific expression is given in the appendix [29].Then, combining equations (22), (31), (32), and (34), C2t can be directly calculated as follows:
(35)C2t=h6μ4TfΦfG−1∫E0EKεdεG−1TΦTfTTf.Finally, according to equations (19), (24), (25), and (35), the spacecraft orbit error distribution at any time can be obtained directly by analytical calculation:
(36)mt=m1t+m2t=06×1,Ct=C1t+C2t=φt,t0C0LIφTt,t0+h6μ4TfΦfG−1∫E0EKεdεG−1TΦTfTTf.


## 3. The Calculation Method of Space Interception Probability

According to the analytical calculation method for the uncertainty propagation of spacecraft orbits given in Section2, 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 TCAFor 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 Section2).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 radiusWhen 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 matrixCptTCA at the TCA will be used in the calculation. Divide CtTCA as shown in Figure 2, then
(37)CptTCA=C11C12C13C21C22C23C31C32C33.Figure 2
Partition ofCt.At the TCA, the relative position vector between the target distribution center and the interceptor distribution center isrTCA. According to the assumptions 1, 2, and 3, it can be obtained that the rTCA and the relative velocity vTCA are perpendicular to each other at the TCA. The plane that is perpendicular to vTCA and passes through rTCA 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 targetO1 as the origin, the xe-axis points to the projection point of the interceptor O2 in the encounter plane, the ze-axis points to the vTCA direction, and the ye-axis satisfies the right-hand rule (as shown in Figure 3). The unit vectors of the three axes can be expressed as follows:
(38)ie=rTCArTCA,ke=vTCAvTCA,je=ke×ie.Figure 3
The encounter coordinate system.Then, the transition matrix from the Earth Centered Inertial (ECI) coordinate system to the encounter coordinate system is(39)MECI→ECS=iejekeT.Let the transition matrix from the LVLHI coordinate system to the ECI coordinate system beMLI→ECI, then the position vectors of the target and the interceptor in the encounter coordinate system are
(40)Re1=MECI→ECSMLI1→ECIR1,Re2=MECI→ECSMLI2→ECIR2,where Re1 and Re2 are the position vectors of the target and the interceptor in the encounter coordinate system, respectively, and R1 and R2 are the position vectors of the target and the interceptor in the LVLHI coordinate system, respectively.In Figure3, O1xeye is the encounter plane between the target and the interceptor.


### 3.2. Interception Probability Calculation Model

Through the calculation of error propagation in Section2, 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:
(41)Cpe1tTCA=MECI→ECSMLI1→ECICp1tTCAMTLI1→ECIMTECI→ECS,Cpe2tTCA=MECI→ECSMLI2→ECICp2tTCAMTLI2→ECIMTECI→ECS,where Cp1tTCA and Cp2tTCA are, respectively, the position error covariance matrices of the target and the interceptor in the LVLHI coordinate system at the TCA and Cpe1tTCA and Cpe2tTCA 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 beRep1 and Rep2, respectively. So, the mean value of the relative position vector rep between the distribution centers of the target and interceptor in the encounter plane can be expressed as follows:
(42)Erep=ERep2−Rep1=retTCA0T,where re 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 beCpep1tTCA and Cpep2tTCA, respectively. Then, refer to the composition of Ct (as shown in Figure 2), the expressions of Cpep1tTCA and Cpep2tTCA are as follows:
(43)Cpep1tTCA=Cpe111Cpe112Cpe121Cpe122,Cpep2tTCA=Cpe211Cpe212Cpe221Cpe222,where Cpe1ij and Cpe2ij represent the elements at the positions of the ith row and the jth column in the Cpe1tTCA and Cpe2tTCA, respectively, Cpe112 is equal to Cpe121, and Cpe212 is equal to Cpe221.Then, at the TCA, the error covariance matrix of the relative position vectorrep between the target and the interceptor in the encounter plane is
(44)CpeptTCA=varrep=varRep2−Rep1=Cpep2tTCA+Cpep1tTCA−2covRep2,Rep1=Cpep2tTCA+Cpep1tTCA.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 d1 and d2, 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 d=d1+d2). 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 Figure4 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), CpeptTCA is a symmetric matrix. So, the joint error covariance matrix in the encounter plane at the TCA is
(45)CpeptTCA=Cxe2κCxeCyeκCxeCyeCye2,where κ∈−1,1 is the correlation coefficient. Cxe,Cye∈ℝ+ are the standard deviations of the joint error ellipse along the xe-axis and ye-axis, respectively.It can be seen from Figure4 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 x-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 O1xcyc. Figure 5 shows the rotation process.Figure 5
Rotation of the encounter coordinate system to the calculation coordinate system.In Figure5, mxc and myc are the expected values of the joint error ellipse along the xc-axis and yc-axis directions, respectively. θ∈−π/2,π/2 is the rotation angle, and its calculation method is discussed in the following three cases:
(1)
Ifκ=0 and Cxe≥Cye, then θ=0
(2)
Ifκ=0 and Cxe<Cye, then θ=π/2
(3)
Ifκ≠0, then θ is calculated by
(46)θ=−arctanCye2−Cxe22κCxeCye+signκ1+Cye2−Cxe22κCxeCye2.After rotation, the joint error covariance matrixCpctTCAin the calculation coordinate system is
(47)CpctTCA=Cxc200Cyc2,where Cxc and Cyc are the standard deviations of the joint error ellipse along with the xc-axis and yc-axis directions, respectively. At this time, the distribution parameters of the joint error in the calculation coordinate system are
(48)mxc=recosθ,myc=resinθ,Cxc2=Cxe2cos2θ−κCxeCyesin2θ+Cye2sin2θ,Cyc2=Cxe2sin2θ+κCxeCyesin2θ+Cye2cos2θ.Then, the 2-D PDF can be expressed as follows:(49)fxc,yc=12πCxcCycexp−12xc−mxc2Cxc2+yc−myc2Cyc2.The interception probabilityPc can be expressed as the integral of the PDF in the joint circle domain:
(50)Pc=∬xc2+yc2≤d2/4fxc,ycdxcdyc.


## 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 TCAFor 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 Section2).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 radiusWhen 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 matrixCptTCA at the TCA will be used in the calculation. Divide CtTCA as shown in Figure 2, then
(37)CptTCA=C11C12C13C21C22C23C31C32C33.Figure 2
Partition ofCt.At the TCA, the relative position vector between the target distribution center and the interceptor distribution center isrTCA. According to the assumptions 1, 2, and 3, it can be obtained that the rTCA and the relative velocity vTCA are perpendicular to each other at the TCA. The plane that is perpendicular to vTCA and passes through rTCA 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 targetO1 as the origin, the xe-axis points to the projection point of the interceptor O2 in the encounter plane, the ze-axis points to the vTCA direction, and the ye-axis satisfies the right-hand rule (as shown in Figure 3). The unit vectors of the three axes can be expressed as follows:
(38)ie=rTCArTCA,ke=vTCAvTCA,je=ke×ie.Figure 3
The encounter coordinate system.Then, the transition matrix from the Earth Centered Inertial (ECI) coordinate system to the encounter coordinate system is(39)MECI→ECS=iejekeT.Let the transition matrix from the LVLHI coordinate system to the ECI coordinate system beMLI→ECI, then the position vectors of the target and the interceptor in the encounter coordinate system are
(40)Re1=MECI→ECSMLI1→ECIR1,Re2=MECI→ECSMLI2→ECIR2,where Re1 and Re2 are the position vectors of the target and the interceptor in the encounter coordinate system, respectively, and R1 and R2 are the position vectors of the target and the interceptor in the LVLHI coordinate system, respectively.In Figure3, O1xeye is the encounter plane between the target and the interceptor.


## 3.2. Interception Probability Calculation Model

Through the calculation of error propagation in Section2, 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:
(41)Cpe1tTCA=MECI→ECSMLI1→ECICp1tTCAMTLI1→ECIMTECI→ECS,Cpe2tTCA=MECI→ECSMLI2→ECICp2tTCAMTLI2→ECIMTECI→ECS,where Cp1tTCA and Cp2tTCA are, respectively, the position error covariance matrices of the target and the interceptor in the LVLHI coordinate system at the TCA and Cpe1tTCA and Cpe2tTCA 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 beRep1 and Rep2, respectively. So, the mean value of the relative position vector rep between the distribution centers of the target and interceptor in the encounter plane can be expressed as follows:
(42)Erep=ERep2−Rep1=retTCA0T,where re 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 beCpep1tTCA and Cpep2tTCA, respectively. Then, refer to the composition of Ct (as shown in Figure 2), the expressions of Cpep1tTCA and Cpep2tTCA are as follows:
(43)Cpep1tTCA=Cpe111Cpe112Cpe121Cpe122,Cpep2tTCA=Cpe211Cpe212Cpe221Cpe222,where Cpe1ij and Cpe2ij represent the elements at the positions of the ith row and the jth column in the Cpe1tTCA and Cpe2tTCA, respectively, Cpe112 is equal to Cpe121, and Cpe212 is equal to Cpe221.Then, at the TCA, the error covariance matrix of the relative position vectorrep between the target and the interceptor in the encounter plane is
(44)CpeptTCA=varrep=varRep2−Rep1=Cpep2tTCA+Cpep1tTCA−2covRep2,Rep1=Cpep2tTCA+Cpep1tTCA.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 d1 and d2, 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 d=d1+d2). 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 Figure4 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), CpeptTCA is a symmetric matrix. So, the joint error covariance matrix in the encounter plane at the TCA is
(45)CpeptTCA=Cxe2κCxeCyeκCxeCyeCye2,where κ∈−1,1 is the correlation coefficient. Cxe,Cye∈ℝ+ are the standard deviations of the joint error ellipse along the xe-axis and ye-axis, respectively.It can be seen from Figure4 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 x-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 O1xcyc. Figure 5 shows the rotation process.Figure 5
Rotation of the encounter coordinate system to the calculation coordinate system.In Figure5, mxc and myc are the expected values of the joint error ellipse along the xc-axis and yc-axis directions, respectively. θ∈−π/2,π/2 is the rotation angle, and its calculation method is discussed in the following three cases:
(1)
Ifκ=0 and Cxe≥Cye, then θ=0
(2)
Ifκ=0 and Cxe<Cye, then θ=π/2
(3)
Ifκ≠0, then θ is calculated by
(46)θ=−arctanCye2−Cxe22κCxeCye+signκ1+Cye2−Cxe22κCxeCye2.After rotation, the joint error covariance matrixCpctTCAin the calculation coordinate system is
(47)CpctTCA=Cxc200Cyc2,where Cxc and Cyc are the standard deviations of the joint error ellipse along with the xc-axis and yc-axis directions, respectively. At this time, the distribution parameters of the joint error in the calculation coordinate system are
(48)mxc=recosθ,myc=resinθ,Cxc2=Cxe2cos2θ−κCxeCyesin2θ+Cye2sin2θ,Cyc2=Cxe2sin2θ+κCxeCyesin2θ+Cye2cos2θ.Then, the 2-D PDF can be expressed as follows:(49)fxc,yc=12πCxcCycexp−12xc−mxc2Cxc2+yc−myc2Cyc2.The interception probabilityPc can be expressed as the integral of the PDF in the joint circle domain:
(50)Pc=∬xc2+yc2≤d2/4fxc,ycdxcdyc.


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

In Section3, 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 Pc=gλ, where λ=d2/4. The form of function g:ℝ+↦ℝ+ is as follows:
(51)gλ=12πCxcCyc∬xc2+yc2≤λ2exp−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdyc.Since each of the terms obtained by the functiongλ 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 q should be set before the operation. Replacing g with q⋅g can avoid the elimination phenomenon. From the perspective of the Laplace transform, we use the form of function q as q=epλ, where p is the regulatory factor. According to the literature [24, 30], the selectable factor p is 1/2Cyc2.According to the Laplace transform theorem, the Laplace transform of functiongλ can be expressed as follows:
(52)Lgs=∫0+∞gλexp−sλdλ=12πCxcCyc⋅∫0+∞∬xc2+yc2≤λ2exp−sλ−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdycdλ=12πCxcCycs∫ℝ2exp−sxc2+yc2−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdyc=12πCxcCycs⋅∫ℝ2exp−sxc2−12xc−mxc2Cxc2exp−syc2−12yc−myc2Cyc2dxcdyc,where s is a complex variable.According to the Fubini theorem [31], when fxc,yc=f1xcf2yc, there is
(53)∫ℝ2fxc,ycdxcdyc=∫−∞+∞f1xcdxc∫−∞+∞f2ycdyc.And the final expression forLgs can be
(54)Lgs=exp−smxc2/2sCxc2+1+myc2/2sCyc2+1s2sCxc2+12sCyc2+1.Define the functionhλ:ℝ+↦ℝ+ as follows:
(55)hλ=epλgλ.Then, the Laplace transform of functionhλ can be expressed as follows:
(56)Lhs=∫0+∞gλexppλexp−sλdλ=∫0+∞gλexp−s−pλdλ=Lgs−p=exp−s−pmxc2/2s−pCxc2+1+myc2/2s−pCyc2+1s−p2s−pCxc2+12s−pCyc2+1,where s>p. Convert equation (56) to
(57)Lhs=b0expQx/s−pη+Qy/sss−pηs−p,where the expressions of η, Qx, Qy, and b0 are as follows:
(58)η=1−Cyc2Cxc2,Qx=mxc24Cxc4,Qy=myc24Cyc4,b0=12CxcCycexp−12mxc2Cxc2+myc2Cyc2.Next, we perform the Taylor expansion onLhs 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 ofLhs−1 be
(59)Lhs−1=s2b0expsQx/1−pηs+sQy1−pηs1−ps.DefineL~hs=s−2Lhs−1, then the derivative of L~hs is
(60)dL~hds=b0expsQx/1−pηs+sQy1−pηs1−ps2Qx1−pηs2+Qy⋅1−pηs1/21−ps−pηps−121−pηs−p1−pηs1/2=L~hsQy+Qx1−pηs2+pη21−pηs+p1−ps.According toL~h0=b0 and equation (60), L~hs is differentiable in the complex plane. According to the definition of power series, L~hs can be expanded into the form of power series:
(61)L~hs=∑n=0∞bnsn,where
(62)b0=12CxcCycexp−12mxc2Cxc2+myc2Cyc2,b1=b0p1+η2+Qx+Qyp,b2=b02p21+η2+Qx+Qyp2+η22+1+2ηQxp,b3=b06p31+η2+Qx+Qyp3+31+η2+Qx+Qyp⋅η22+1+2ηQxp+2η32+1+3η2Qxp,bn+4=−p3η2Qybnn+4+p2ηpηn+2.5+2Qyη/2+1bn+1n+4−ppηη/2+12n+5+η2Qy+1.5p+Qx+Qybn+2n+4+p2η+1n+3+pη/2+1+Qx+Qybn+3n+4.SinceL~hs=s−2Lhs−1, the power series expression of Lhs is
(63)Lhs=1s2∑n=0∞bn1sn=∑n=0∞bnsn+2.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:
(64)Pc=gλ=e−pλhλ=e−pλLh−1s=e−pd2/4∑n=0+∞bnd2n+24n+1n+1!.Letcn=bnd2n+2/4n+1n+1!, then (63) can be written as follows:
(65)Pc=e−pd2/4∑n=0+∞cn,where
(66)c0=b0d24,c1=b0d4p321+η2+Qx+Qyp,c2=b0d6p27681+η2+Qx+Qyp2+η22+1+2ηQxp,c3=b0d8p3368641+η2+Qx+Qyp3+31+η2+Qx+Qyp⋅η22+1+2ηQxp+2η32+1+3η2Qxp,cn+4=−d8p3η2Qycn256n+2n+3n+42n+5+d6p2ηpηn+2.5+2Qyη/2+1cn+164n+3n+42n+5−d4ppηη/2+12n+5+η2Qy+1.5p+Qx+Qycn+216n+42n+5+d2p2η+1n+3+pη/2+1+Qx+Qycn+34n+4n+5.If the sum of the firstn 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 n terms of the power series be Pc′, then the truncation error Sn is
(67)Sn=Pc−Pc′=e−pd2/4∑ℓ=n+∞bℓd2ℓ+24ℓ+1ℓ+1!.According to equation (67), when the number n of terms of the power series and the diameter d of the joint domain is given, Sn is only related to bℓ. Define the variables b¯ℓ and b~ℓ as follows:
(68)b¯ℓ=b0pℓ1+η2+Qx+Qypℓ,b~ℓ=b0pℓ.According to equation (62), the following relationship can be obtained:
(69)n+1bn+1=∑ℓ=0nχℓbn−ℓ,where the expression of χℓ is as follows:
(70)χℓ=pℓ+11+ηℓQxℓ+1p+η2+0ℓ>0Qyℓ=0.Sincep, η, Qx, and Qy are nonnegative, we can get
(71)b0χℓb~ℓ+1≥1≥b0χℓb¯ℓ+1.According to equation (71), we can get
(72)n+1b~n+1≤b~n+1∑ℓ=0nb0χℓb~ℓ+1=∑ℓ=0nχℓb~n−ℓ,n+1b¯n+1≥b¯n+1∑ℓ=0nb0χℓb¯ℓ+1=∑ℓ=0nχℓb¯n−ℓ.Therefore,b~ℓ≤bℓ≤b¯ℓ can be easily obtained by proof by contradiction and the value range of truncation error Sn in the calculation of space interception probability is
(73)e−pd2/4∑ℓ=n+∞b~ℓd2ℓ+24ℓ+1ℓ+1!≤Sn≤e−pd2/4∑ℓ=n+∞b¯ℓd2ℓ+24ℓ+1ℓ+1!.If℘>0, then the following relations hold
(74)℘d2/4n+1℘n+1!≤∑ℓ=n∞℘ℓd2/4ℓ+1ℓ+1!≤e℘d2/4℘d2/4n+1℘n+1!.Therefore, equation (73) can be rewritten as follows:
(75)Sn≥b0e−pd2/4pn+1d2n+24n+1pn+1!,Sn≤b0exppd2/4η/2+Qx+Qy/pp1+η/2+Qx+Qy/pn+1!⋅pd241+η2+Qx+Qypn+1=S¯,where S¯ is the maximum truncation error.Therefore, by combining Sections2 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.
1. Input:
Orbital information and orbital error distribution of target and interceptor at the initial momentXTLIt0∼ℕXTLIt0;XT0LI,CT0LI, XILIt0∼ℕXILIt0;XI0LI,CI0LI; Target and interceptor diameter d1, d2; 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 anomalyfT of the target orbit and the true anomaly fI of the interceptor orbit corresponding to the TCA;
3) CalculateCT1tTCA and CI1tTCA according to formula (25);
4) Calculate the eccentric anomalyET of the target orbit and the eccentric anomaly EI of the interceptor orbit corresponding to the TCA;
5) CalculateCT2tTCA, CI2tTCA according to the appendix, formula (6), formula (30), and formula (35);
6) According to equation (36), the state error covariance matrix CTtTCA and CItTCA of the target and the interceptor at the TCA are obtained from CT1tTCA+CT2tTCA and CI1tTCA+CI2tTCA.
3. Error projection and coordinate system rotation
1) According to equations (41), (43), and (44), calculate the joint error covariance matrix CpeptTCA 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 mxc, myc, Cxc, and Cyc in the calculation coordinate system according to equation (48).
4. Space interception probability calculation
1) Calculateη, Qx, Qy, and b0 according to equation (58), and calculate c0~cN according to equation (66) and the term number N of power series;
2) Calculate space interception probabilityPc and maximum truncation error S¯ according to equations (65) and (75).
5. Output:
Space interception probabilityPc and maximum truncation error S¯.


## 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 Pc=gλ, where λ=d2/4. The form of function g:ℝ+↦ℝ+ is as follows:
(51)gλ=12πCxcCyc∬xc2+yc2≤λ2exp−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdyc.Since each of the terms obtained by the functiongλ 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 q should be set before the operation. Replacing g with q⋅g can avoid the elimination phenomenon. From the perspective of the Laplace transform, we use the form of function q as q=epλ, where p is the regulatory factor. According to the literature [24, 30], the selectable factor p is 1/2Cyc2.According to the Laplace transform theorem, the Laplace transform of functiongλ can be expressed as follows:
(52)Lgs=∫0+∞gλexp−sλdλ=12πCxcCyc⋅∫0+∞∬xc2+yc2≤λ2exp−sλ−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdycdλ=12πCxcCycs∫ℝ2exp−sxc2+yc2−12xc−mxc2Cxc2+yc−myc2Cyc2dxcdyc=12πCxcCycs⋅∫ℝ2exp−sxc2−12xc−mxc2Cxc2exp−syc2−12yc−myc2Cyc2dxcdyc,where s is a complex variable.According to the Fubini theorem [31], when fxc,yc=f1xcf2yc, there is
(53)∫ℝ2fxc,ycdxcdyc=∫−∞+∞f1xcdxc∫−∞+∞f2ycdyc.And the final expression forLgs can be
(54)Lgs=exp−smxc2/2sCxc2+1+myc2/2sCyc2+1s2sCxc2+12sCyc2+1.Define the functionhλ:ℝ+↦ℝ+ as follows:
(55)hλ=epλgλ.Then, the Laplace transform of functionhλ can be expressed as follows:
(56)Lhs=∫0+∞gλexppλexp−sλdλ=∫0+∞gλexp−s−pλdλ=Lgs−p=exp−s−pmxc2/2s−pCxc2+1+myc2/2s−pCyc2+1s−p2s−pCxc2+12s−pCyc2+1,where s>p. Convert equation (56) to
(57)Lhs=b0expQx/s−pη+Qy/sss−pηs−p,where the expressions of η, Qx, Qy, and b0 are as follows:
(58)η=1−Cyc2Cxc2,Qx=mxc24Cxc4,Qy=myc24Cyc4,b0=12CxcCycexp−12mxc2Cxc2+myc2Cyc2.Next, we perform the Taylor expansion onLhs 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 ofLhs−1 be
(59)Lhs−1=s2b0expsQx/1−pηs+sQy1−pηs1−ps.DefineL~hs=s−2Lhs−1, then the derivative of L~hs is
(60)dL~hds=b0expsQx/1−pηs+sQy1−pηs1−ps2Qx1−pηs2+Qy⋅1−pηs1/21−ps−pηps−121−pηs−p1−pηs1/2=L~hsQy+Qx1−pηs2+pη21−pηs+p1−ps.According toL~h0=b0 and equation (60), L~hs is differentiable in the complex plane. According to the definition of power series, L~hs can be expanded into the form of power series:
(61)L~hs=∑n=0∞bnsn,where
(62)b0=12CxcCycexp−12mxc2Cxc2+myc2Cyc2,b1=b0p1+η2+Qx+Qyp,b2=b02p21+η2+Qx+Qyp2+η22+1+2ηQxp,b3=b06p31+η2+Qx+Qyp3+31+η2+Qx+Qyp⋅η22+1+2ηQxp+2η32+1+3η2Qxp,bn+4=−p3η2Qybnn+4+p2ηpηn+2.5+2Qyη/2+1bn+1n+4−ppηη/2+12n+5+η2Qy+1.5p+Qx+Qybn+2n+4+p2η+1n+3+pη/2+1+Qx+Qybn+3n+4.SinceL~hs=s−2Lhs−1, the power series expression of Lhs is
(63)Lhs=1s2∑n=0∞bn1sn=∑n=0∞bnsn+2.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:
(64)Pc=gλ=e−pλhλ=e−pλLh−1s=e−pd2/4∑n=0+∞bnd2n+24n+1n+1!.Letcn=bnd2n+2/4n+1n+1!, then (63) can be written as follows:
(65)Pc=e−pd2/4∑n=0+∞cn,where
(66)c0=b0d24,c1=b0d4p321+η2+Qx+Qyp,c2=b0d6p27681+η2+Qx+Qyp2+η22+1+2ηQxp,c3=b0d8p3368641+η2+Qx+Qyp3+31+η2+Qx+Qyp⋅η22+1+2ηQxp+2η32+1+3η2Qxp,cn+4=−d8p3η2Qycn256n+2n+3n+42n+5+d6p2ηpηn+2.5+2Qyη/2+1cn+164n+3n+42n+5−d4ppηη/2+12n+5+η2Qy+1.5p+Qx+Qycn+216n+42n+5+d2p2η+1n+3+pη/2+1+Qx+Qycn+34n+4n+5.If the sum of the firstn 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 n terms of the power series be Pc′, then the truncation error Sn is
(67)Sn=Pc−Pc′=e−pd2/4∑ℓ=n+∞bℓd2ℓ+24ℓ+1ℓ+1!.According to equation (67), when the number n of terms of the power series and the diameter d of the joint domain is given, Sn is only related to bℓ. Define the variables b¯ℓ and b~ℓ as follows:
(68)b¯ℓ=b0pℓ1+η2+Qx+Qypℓ,b~ℓ=b0pℓ.According to equation (62), the following relationship can be obtained:
(69)n+1bn+1=∑ℓ=0nχℓbn−ℓ,where the expression of χℓ is as follows:
(70)χℓ=pℓ+11+ηℓQxℓ+1p+η2+0ℓ>0Qyℓ=0.Sincep, η, Qx, and Qy are nonnegative, we can get
(71)b0χℓb~ℓ+1≥1≥b0χℓb¯ℓ+1.According to equation (71), we can get
(72)n+1b~n+1≤b~n+1∑ℓ=0nb0χℓb~ℓ+1=∑ℓ=0nχℓb~n−ℓ,n+1b¯n+1≥b¯n+1∑ℓ=0nb0χℓb¯ℓ+1=∑ℓ=0nχℓb¯n−ℓ.Therefore,b~ℓ≤bℓ≤b¯ℓ can be easily obtained by proof by contradiction and the value range of truncation error Sn in the calculation of space interception probability is
(73)e−pd2/4∑ℓ=n+∞b~ℓd2ℓ+24ℓ+1ℓ+1!≤Sn≤e−pd2/4∑ℓ=n+∞b¯ℓd2ℓ+24ℓ+1ℓ+1!.If℘>0, then the following relations hold
(74)℘d2/4n+1℘n+1!≤∑ℓ=n∞℘ℓd2/4ℓ+1ℓ+1!≤e℘d2/4℘d2/4n+1℘n+1!.Therefore, equation (73) can be rewritten as follows:
(75)Sn≥b0e−pd2/4pn+1d2n+24n+1pn+1!,Sn≤b0exppd2/4η/2+Qx+Qy/pp1+η/2+Qx+Qy/pn+1!⋅pd241+η2+Qx+Qypn+1=S¯,where S¯ is the maximum truncation error.Therefore, by combining Sections2 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.
1. Input:
Orbital information and orbital error distribution of target and interceptor at the initial momentXTLIt0∼ℕXTLIt0;XT0LI,CT0LI, XILIt0∼ℕXILIt0;XI0LI,CI0LI; Target and interceptor diameter d1, d2; 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 anomalyfT of the target orbit and the true anomaly fI of the interceptor orbit corresponding to the TCA;
3) CalculateCT1tTCA and CI1tTCA according to formula (25);
4) Calculate the eccentric anomalyET of the target orbit and the eccentric anomaly EI of the interceptor orbit corresponding to the TCA;
5) CalculateCT2tTCA, CI2tTCA according to the appendix, formula (6), formula (30), and formula (35);
6) According to equation (36), the state error covariance matrix CTtTCA and CItTCA of the target and the interceptor at the TCA are obtained from CT1tTCA+CT2tTCA and CI1tTCA+CI2tTCA.
3. Error projection and coordinate system rotation
1) According to equations (41), (43), and (44), calculate the joint error covariance matrix CpeptTCA 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 mxc, myc, Cxc, and Cyc in the calculation coordinate system according to equation (48).
4. Space interception probability calculation
1) Calculateη, Qx, Qy, and b0 according to equation (58), and calculate c0~cN according to equation (66) and the term number N of power series;
2) Calculate space interception probabilityPc and maximum truncation error S¯ according to equations (65) and (75).
5. Output:
Space interception probabilityPc and maximum truncation error S¯.


## 5. Simulation Example

In this section, a simulation example is given to illustrate the correctness of the proposed algorithm. Take the gravitational constantμ for 398600.4418km3/s2, 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 axisakm
16545
36834.9
Eccentricitye
0.4
0.676229
Orbit inclinationi°
63
114.433
RAANΩ°
25
191.606
Argument of perigeeϖ°
20
111.057
Initial true anomalyf°
30
287.59Select 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). Figure6 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. Tables2 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
x (km)
y (km)
z (km)
Target
-4291.771159
4390.710760
11366.515836
Interceptor
-4291.821845
4390.691209
11366.550420Table 3
Velocity parameters of the target and interceptor at the TCA.
Velocity in the J2000 coordinate system
vx (km/s)
vy (km/s)
vz (km/s)
Target
-5.681314
-2.060952
1.042846
Interceptor
6.981309
1.460946
0.057147At 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 Section2, the comparative analysis is performed by Monte Carlo simulation. Suppose the initial error distribution of the target and the interceptor is
(76)XT0LI=06×6XI0LI=06×6,CT0LI=diag302m2302m2302m20.012m2/s20.012m2/s20.012m2/s2,CI0LI=diag102m2102m2102m20.12m2/s20.12m2/s20.12m2/s2.According to the literature [32], the intensity of U is
(77)U=diag10−510−510−5m2/s4.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 Figures7 and 8.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 thex-axis direction in the LVLHI coordinate system. (b) Position error mean square error propagation curves of the target and interceptor in the y-axis direction in the LVLHI coordinate system. (c) Position error mean square error propagation curves of the target and interceptor in the z-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 thex-axis direction in the LVLHI coordinate system. (b) Velocity error mean square error propagation curves of the target and interceptor in the y-axis direction in the LVLHI coordinate system. (c) Velocity error mean square error propagation curves of the target and interceptor in the z-axis direction in the LVLHI coordinate system.
(a)
(b)
(c)It can be seen from Figures7 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%. Figure10 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 Figures9 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.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.
(a)
(b)
(c)From Figures11 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
(78)CpeptTCA=0.02519−0.00325−0.003250.02362km2,of the target and the interceptor in the encounter plane at the TCA can be calculated. Then, according to equation (46), the rotation angle θ=38.21057° of the coordinate system can be calculated. Finally, the joint error distribution parameters mxc=0.02427km, myc=0.01911km, Cxc2=0.02775km2, and Cyc2=0.02106km2 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 η=0.24123, Qx=0.19127km−2, Qy=0.20589km−2, and b0=20.28845km−2can be calculated. Take the power series terms number N=1,2,⋯,8. According to equations (65) and (75), the space interception probability Pc and maximum truncation error S¯ are calculated, respectively. Table 4 shows the calculation results.Table 4
Space interception probabilityPc and maximum truncation error S¯ through the Laplace transform.
N
Pc
S¯
1
5.273768358×10−4
3.089214758×10−8
2
5.275620585×10−4
1.084979020×10−11
3
5.275621365×10−4
2.540407118×10−15
4
5.275621366×10−4
4.461147316×10−19
5
5.275621366×10−4
6.267290068×10−23
6
5.275621366×10−4
7.337224040×10−27
7
5.275621366×10−4
7.362697802×10−31
8
5.275621366×10−4
6.464727505×10−35It can be seen from Table4 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 N 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 N 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 probabilityPc obtained by the three methods.
Pc
Monte Carlo
5.274333333×10−4
Chan
5.275634633×10−4
Laplace transform
5.275621366×10−4It can be seen from Table5 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.67015Table6 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. Figure14 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 is5.3295×10−4. 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



---

*Source: 5316396-2019-11-03.xml*