Portal de Periódicos da CAPES

Satellite Instantaneous Collision Probability Computation Using Equivalent Volume Cuboids

2020; American Institute of Aeronautics and Astronautics; Volume: 43; Issue: 9 Linguagem: Inglês

10.2514/1.g004711

1533-3884

Shuo Zhang, Tuo Fu, Defeng Chen, Huawei Cao,

Risk and Safety Analysis

Open AccessEngineering NotesSatellite Instantaneous Collision Probability Computation Using Equivalent Volume CuboidsShuo Zhang, Tuo Fu, Defeng Chen and Huawei CaoShuo ZhangBeijing Institute of Technology, 100081 Beijing, People's Republic of China, Tuo FuBeijing Institute of Technology, 100081 Beijing, People's Republic of China, Defeng ChenBeijing Institute of Technology, 100081 Beijing, People's Republic of China and Huawei CaoBeijing Institute of Technology, 100081 Beijing, People's Republic of ChinaPublished Online:21 Jun 2020https://doi.org/10.2514/1.G004711SectionsRead Now ToolsAdd to favoritesDownload citationTrack citations ShareShare onFacebookTwitterLinked InRedditEmail AboutI. IntroductionSatellite conjunction analysis is a fundamental task for space situational awareness (SSA). Because of the rapid growth of launch activities over the past several decades, space is becoming increasingly congested. As a consequence, the risk of collisions among resident space objects (RSOs) has also increased. To protect high-value satellites or the manned space station, an operator must determine whether there is a need for a maneuver to avoid any underlying collisions. Because of the imperfect knowledge of target positions and velocities, collision risks are usually evaluated in a probabilistic manner. The satellite collision probability is a quantity that can provide the operator with the necessary information to assess the risk level. Therefore, an accurate and effective algorithm to compute the satellite collision probability is of practical value.The problem of accurate conjunction analysis or collision probability evaluation has been extensively investigated. In general, the collision probability can be divided into two categories: the cumulative [1] and the instantaneous collision probabilities. The former is the overall collision probability for two targets during a period of time of interest, whereas the latter refers to the collision probability at a single time instant.Great efforts have been made to quantify the cumulative collision probability in the SSA community, including works by Foster and Estes [2], Chan [3], Alfano [4], Patera [5–7], and the more recent work by Serra et al. [8] and García-Pelayo and Hernando-Ayuso [9]. These methods were all designed for short-term encounters with a Gaussian uncertainty, which led to a two-dimensional (2-D) Gaussian probability density function (PDF) integral over a circle. Among the above methods, Foster's, Alfano's, and Patera's rely on numerical integrations, whereas Chan's, Serra's, and García-Pelayo's are analytic. Moreover, Chan's method is approximate, whereas Serra's and García-Pelayo's methods are exact. However, none of these works considered the velocity uncertainty. Coppola [10] developed an approach to incorporate the velocity uncertainty, which is suitable for both the short-term and long-term encounters. It was shown that the classical short-term encounter result is a special case when there is no uncertainty in the velocity. Because of the nonlinearity of orbital dynamics, the target state with an initial Gaussian distributed uncertainty will become non-Gaussian after the long-term propagation. DeMars et al. [11] used a Gaussian mixture model (GMM) to propagate the orbit uncertainty, which resulted in a GMM representation of the relative state vector of the involved target pair. A numerical integration approach was then adopted to obtain the collision probability. Vittaldev and Russell [12] also adopted the GMM as a tool to quantify the uncertainty at the closest approach time, and the all-on-all evaluation of the collision probability between the two GMMs was conducted both with the linear collision method and Coppola's method. Armellin et al. [13] and Morselli et al. [14,15] used the differential algebra (DA) technique to expand the closest approach distance of two targets as a high-order Taylor series with respect to their initial conditions, and a Monte Carlo simulation was then used to obtain the collision probability. The DA technique transforms the time-consuming numerical integration into a fast polynomial evaluation, leading to a significant reduction in the computation time.Compared with the cumulative collision probability, much less attention has been given to the problems of the instantaneous collision probability. In [16], Chan investigated this problem under the Gaussian uncertainty assumption. It was demonstrated that the core of the instantaneous collision probability computation is the integral of a general three-dimensional (3-D) Gaussian PDF over a sphere. The author subsequently proposed two methods to approximate such an integral. One method was a semi-analytical approach based on the equivalent area, and the other was a fully analytical approach based on the equivalent volume. However, the former required numerical integration along a certain dimension, which could be time-consuming, and the latter was proven to be inaccurate. Furthermore, the Gaussian assumption could be questionable because the state PDF of an RSO would be non-Gaussian after propagation. To adapt to the non-Gaussian distribution of the target state, Jones et al. [17] and Jones and Doostan [18] combined the polynomial chaos expansion (PCE) technique with a Monte Carlo simulation to quantify the satellite instantaneous collision probability. For each target, the authors first used the PCE to construct a surrogate model of the target state at a future time, which was a polynomial of random inputs, and then used a Monte Carlo simulation to estimate the collision probability. However, the Monte Carlo simulation could be inefficient, especially when the collision probability was relatively small. In addition, Adurthi and Singla [19] used the recently proposed conjugate unscented transform (CUT) to estimate the statistical moments of the relative position vector. Its PDF was then approximately reconstructed with these moments using the maximum entropy method. Finally the time-consuming numerical integration was adopted to obtain the collision probability.In this paper, we propose an integrated, two-stage approach to effectively compute the satellite instantaneous collision probability. First, it uses the GMM to approximate the relative position vector PDF for the two involved RSOs at a future time, and the GMM is then integrated over the hard body sphere. The core of this integration involves evaluating the integral for a general 3-D Gaussian PDF over a sphere, in which we adopt the equivalent volume method to yield an analytic approximation. The key difference from the adjoining parallelepipeds method in [20] is that we use a single cuboid rather than a bundle of parallelepipeds to approximate the transformed hard body shape. To reduce the approximation error, we force the volume of the cuboid to equal the volume of the transformed hard body ellipsoid. Using a single cuboid is beneficial for reducing the computation cost. The main contributions of our work are as follows: In the first stage, we propose a new split strategy to effectively represent the initial target Gaussian PDF as a GMM. We select the most significant eigendirection each time to apply the splitting library rather than the conventional greatest nonlinear direction [12] or the spectral direction corresponding to the maximum eigenvalue [21] of the covariance matrix. The most significant eigendirection is along one principle axis of the error ellipsoid, and it is determined based on the linearized position uncertainty analysis.In the second stage, the general Gaussian PDF is mapped into the standard Gaussian PDF, whereas the hard body sphere is mapped into an ellipsoid after applying a series of affine transforms. This ellipsoid is then replaced with an equivalent volume cuboid rather than the equivalent volume sphere. If the ellipsoid is elongated or disk-shaped, the overlapped volume of the cuboid with this ellipsoid will be larger than that of the sphere, leading to an improvement in the accuracy. For a space object, because its position uncertainty along the in-track [22] direction grows faster than along the other directions, it is often the case that the transformed ellipsoid is elongated or disk-like. Hence, our method is expected to perform better than the equivalent volume sphere method in most cases. Furthermore, the equivalent volume cuboid approximation fully decouples the 3-D multiple integration into the product of three one-dimensional (1-D) integrations with explicit expressions. Therefore, there is no need for a second approximation to yield the final integration result.The rest of this paper is organized as follows. Section II presents the definition of the satellite instantaneous collision probability. Section III deals with the uncertainty propagation problem using the GMM approach. Section IV focuses on the equivalent volume cuboid approximation for the GMM integral over the hard body sphere. Section V demonstrates the numerical simulation results for two distinct scenarios. Finally, Sec. VI draws the conclusions of this paper.II. Definition of Satellite Instantaneous Collision ProbabilityThe basics for SSA are various types of observations as obtained from sensors such as radars and optical facilities. These observations are inevitably corrupted by random noise. As a result, the state vector of a space target, which is estimated from noisy observations, is also a random vector. At a given time t, suppose that the relative position vector of two space targets is s with a PDF of p(s). The instantaneous collision probability Pc(t) is then defined as Pc(t)=∭Vp(s) ds (1)where the integration volume V (also referred to as the hard body shape) is fully determined by the shapes of the two targets. If the two targets are modeled as spheres with radii r1 and r2, respectively, then the hard body shape V is also a sphere with radius r=r1+r2. The reason for modeling the two targets as spheres is that it does not require precise knowledge for the shapes and attitudes of the two targets, which greatly simplifies the analysis and simulation processes. Under such an assumption, only the PDF p(s) evolves with time.III. Uncertainty Propagation Using the Gaussian Mixture ModelA central issue for the accurate conjunction analysis is to propagate the initial orbit uncertainty to the desired time t. The time durations between the epoch t0 and time instants with high collision probabilities can be relatively long, which typically vary from several hours to several days. However, under such long periods, the state vector PDF can become non-Gaussian due to the nonlinearity of the orbital dynamics.The GMM is a powerful tool to represent the PDF for the state of a dynamical system at any given time t. It uses a series of weighted Gaussian PDFs to approximate the true PDF as follows: p(x)≈∑i=1MwiN(x;μi,Pi)(2)where x is the state vector of the dynamical system and M is the number of components in the GMM. Each Gaussian component can be fully characterized as a triplet (wi,μi,Pi), where wi, μi, and Pi are the weight, mean vector, and covariance matrix, respectively.It has been proven that a GMM can converge uniformly to the true PDF when the number of components M tends to infinity and the covariance matrix of each component Pi tends to the zero matrix [23].In this section, we propose an integrated method to propagate the orbital uncertainty. The main contribution of our method is the selection of the splitting direction. We select the most significant eigendirection each time rather than the conventional greatest nonlinear direction.A. Split of the Initial Gaussian PDFThe basic question is how to effectively and accurately split the initial state PDF. Although the GMM can converge to the concerned PDF when the number of components becomes infinitely large and each component covariance matrix becomes infinitely small, a good splitting method is indispensable due to the limited computational resources in practical usage.The building block for splitting a multivariate Gaussian PDF is the split of the standard univariate normal PDF N(x;0,1). A commonly used approach [22,23] is based on the optimization method to construct a splitting library that minimizes the Lp distance between N(x;0,1) and the GMM approximation under some constraints. Here, we adopt the method proposed by Horwood and Poore [21], which is based on the Gauss–Hermite quadrature rule. One advantage of such a splitting library is that the statistical moments up to and including the order p for the GMM are equal to those of the original distribution N(0,1).For the sake of simplicity, the initial state PDF of the target is modeled as N(x;μ0,P0), although the GMM method can handle non-Gaussian distributions. A key issue is the selection of the splitting direction. It is noted that a split should not be applied along an arbitrary direction because different directions may have different degrees of nonlinearity. An improper choice of the splitting direction can lead to very limited improvements in the representation accuracy or can even be detrimental in some cases [24,25]. The natural and conventional idea is to select the direction that has the greatest nonlinearity. However, this method requires constructing a function to measure the degrees of nonlinearity along different directions. Here, we present a new approach based on a linear covariance analysis. The linear uncertainty propagator approximates the propagated state as a Gaussian random vector, whose mean vector and covariance matrix are μ=ϕ(t;μ0,t0)(3)P=Φ(t,t0)P0Φ(t,t0)T(4)where ϕ(⋅) is the operator for the numerical integration and Φ(t,t0) is the state transition matrix. Φ(t,t0) can be partitioned as Φ(t,t0)=[ΦrrΦrvΦvrΦvv](5)Then, the linearized position covariance matrix is Prr=[ΦrrΦrv]P0[ΦrrΦrv]T(6)Suppose that P0 can be decomposed as P0=QΛQT(7)where Q is an orthogonal matrix and Λ=diag{λ1,…,λ6} is a diagonal matrix whose elements contain all the eigenvalues of P0. Set Λκ=diag{0,…,0,λκ,0,…,0} and Pκ=QΛκQT. The matrix Pκ can be explained as the covariance matrix of a pseudorandom vector whose uncertainty is solely along the κth symmetrical axis of the error hyperellipsoid defined by P0. Then, the linearly propagated position uncertainty caused by this pseudorandom vector, measured by the position dilution of precision (PDOP), is expressed as σκ=tr{[ΦrrΦrv]Pκ[ΦrrΦrv]T}(8)where tr{⋅} is the trace operator.The splitting direction d is selected as the most significant eigendirection, given by d=maxκ{σκ},κ=1,2,…,6(9)Finally, the unidirectional GMM for the initial state PDF N(x;μ0,P0) is N(x;μ0,P0)≈∑k=1LwkN(x;μ0+QΛμkd,QΛPkdΛTQT)(10)where L is the component number of the splitting library, Λ=diag{λ1,…,λ6}. The elements of μkd are all zeros except for the dth element, which is equal to the kth value in the univariate splitting library. Namely, μkd=[0,…,0,μk,0,…,0]T(11)Pkd is a diagonal matrix whose diagonal elements are all ones except for the dth element, which equals σ2: Pkd=diag{1,…,1,σ2,1,…,1}(12)Equation (10) may not be sufficient to accurately represent the state PDF after propagation, especially when the initial uncertainty is too large, or the propagation duration is too long. In these cases, each component N(x;μ0+QΛμkd,QΛPkdΛTQT) in Eq. (10) can be split similarly as N(x;μ0,P0) to form a multidirectional GMM.The components in the multidirectional GMM grow exponentially with the splitting number. For each component, its weight equals the product of the multiple single-split weights. Most of the GMM components have very small weights after multiple splits because of the unequal weights of the Gauss–Hermite quadrature points. Thus, they have very limited contributions to the final multidirectional GMM representation. To reduce the computational burden, some of the components can be merged using the K-Means clustering algorithm to restrict the final terms for propagation [26].B. Propagation of a Single Gaussian ComponentFor each initial Gaussian term, it is assumed that the Gaussian shape is maintained during the propagation interval [t0,t]. Although there are reported strategies to update the weight to obtain a better representation of the propagated PDF [27,28], we assume the weight to remain unchanged before and after the propagation. Therefore, only the mean vector and the covariance matrix require updating.The linearized propagator expressed by Eqs. (3) and (4) has a superior computation efficiency; however, the propagated covariance matrix could become singular when the propagation duration is long [29]. The scaled unscented transform (SUT) is adopted here to approximate the mean vector and the covariance matrix, which is a variation of the unscented transform (UT) method first proposed by Julier and Uhlmann to overcome the deficiencies of linearization [30–32]. It uses a set of deterministically chosen, weighted sigma points to capture the initial probability distribution.Assuming that the state dimension is N, the process noise can be neglected. For the ith Gaussian component, the corresponding sigma points {ξi,l}l=12N+1, the weights {Wi,lm}l=12N+1, and {Wi,lc}l=12N+1 for mean and covariance are first generated [33], respectively. These initial sigma points are then integrated from t0 to t to yield the propagated sigma points: ηi,l=ϕ(ξi,l,t0,t),l=1,2,…,2N+1(13)Finally, the propagated mean vector and covariance matrix are approximated as follows: μi=∑l=12N+1Wi,lmηi,l(14)Pi=∑l=12N+1Wi,lc(ηi,l−μi)(ηi,l−μi)T(15)C. Verification of the Uncertainty ConsistencyIn practical applications, it is always desirable to determine whether a GMM is an accurate representation of the true PDF. As the most faithful method, Monte Carlo simulation is considered as a good estimation of the true distribution, and the propagated GMM should be compared with Monte Carlo samples. The higher-order Mahalanobis metrics [21] can be used to measure the similarity of the two distributions.First the mean μ1, covariance P1, skewness S1, and fourth-order central moment K1 are estimated using a set of Monte Carlo particles, and μ2, P2, S2, and K2 for the GMM are also obtained [21]. The higher-order Mahalanobis metrics are then calculated as follows: m12=GijΔμiΔμj(16)m22=GikGjlΔPijΔPkl(17)m32=GilGjmGknΔSijkΔSlmn(18)m42=GimGjnGkpGlqΔKijklΔKmnpq(19)where G=(P1+P2)−1, Δμ=μ1−μ2, ΔP=P1−P2, etc. Note that the summation convention is assumed in Eqs. (16–19). All of these metrics are nonnegative and are invariants under affine transforms. Each of them vanishes if the two distributions being compared are the same.Finally, the uncertainty consistency measure (UCM) is defined as UCM=−log(m12+m22+m32+m42)(20)Since mi2,i=1,2,3,4, are invariants under affine transforms, the UCM is also an invariant. A larger UCM value indicates the better consistency. Two identical distributions will result in an infinite UCM value.D. Calculation of the Relative Position Vector PDFThe instantaneous collision probability relies on the relative position vector of the two targets, which can be determined from their marginal position PDFs.Supposing that the initial probabilistic distributions of the two targets are propagated to the same time instant t, let r1 and r2 denote the position vectors for the primary and secondary targets, respectively, and the GMM approximations for r1 and r2 can be easily extracted from their state GMMs as p(r1)≈∑i=1M1w1,iN(r1;μ1,ir,P1,irr)(21)p(r2)≈∑j=1M2w2,jN(r2;μ2,jr,P2,jrr)(22)where μ1,ir and P1,irr are the positional components of the state mean and covariance of the ith GMM term for the first target, respectively. The meanings for μ2,ir and P2,irr can be similarly acquired. It is noted that the right-hand sides of Eqs. (21) and (22) are the exact marginal distributions of their corresponding full-state GMMs, and the Gaussian terms for the two GMMs can be different.The relative position vector s is subsequently computed as s=r1−r2(23)If the initial state vectors of the two targets are mutually independent, the propagated state vectors are also independent and the PDF of s can be calculated as p(s)≈∑i=1M1∑j=1M2w1,iw2,jN(s;μ1,ir−μ2,jr,P1,irr+P2,jrr)(24)IV. Integration of the GMM over a SphereAs described in Sec. II, the computation of the instantaneous collision probability involves integrating the relative position vector PDF over a hard body shape, which is assumed to be a sphere. Because of the additive property of the definite integral, the integration of a GMM over the hard body sphere equals the sum of the integration for each constituent Gaussian PDF over that sphere. Namely, Pc(t)≈∑i=1M1∑j=1M2w1,iw2,j∭VN(s;μ1,ir−μ2,jr,P1,irr+P2,jrr) ds(25)Therefore, the core of Eq. (25) is the evaluation of the integral for a general 3-D Gaussian PDF over a sphere.This section proposes a new analytical method to approximate such an integral. After mapping the general 3-D Gaussian PDF into a standard 3-D Gaussian PDF and the spherical hard body shape into an ellipsoid via a series of affine transforms, we use an equivalent volume cuboid to replace the integration region to yield an analytical approximation for the integral. This equivalent volume cuboid can be viewed as a specialization of Alfano's adjoining parallelepipeds in [20], where the author used a bundle of abutting parallelepipeds to approximate the cumulative collision probability. Figure 1 shows the 2-D illustration of this complete process.Fig. 1 Illustration of the approximation for a Gaussian PDF integral.Consider an arbitrary component in Eq. (25). For simplicity but without loss of generality, let μ and P denote the component mean vector and covariance matrix, respectively. Therefore, the relative position vector s is distributed as N(μ,P). We establish an auxiliary Cartesian coordinate system whose origin is located at the center of the secondary object and whose three axes are aligned with the Earth-centered inertial (ECI) frame. In this coordinate system, the uncertainty can be completely attributed to the secondary object. Denote r as the position vector of the secondary object, and its PDF can be expressed as p(r)=1(2π)3/2|P|1/2exp(−12rTP−1r)(26)The primary object, which is centered at μ, does not contain any uncertainty in its position because we attribute all the uncertainty to the secondary object. The integration region V (the hard body shape) is therefore a sphere centered at μ with radius r, which can be expressed as (r−μ)T(r−μ)≤r2(27)A collision occurs when the secondary object directly enters the hard body sphere. Therefore, in the auxiliary system for the considered component, its collision probability integration is given by I=∭Vp(r) dr(28)This integral is depicted in Fig. 1a, in which the dotted black ellipse and the solid black circle represent the uncertainty ellipse and the hard body sphere, respectively. To our knowledge, there is no exact analytical result reported to solve Eq. (28). Hence we approximate the integration using the following steps.1. Contraction TransformThe covariance matrix is decomposed as P=LLT, and then let r=Lr′. Under this transformation, the original PDF becomes a standard 3-D Gaussian PDF: p(r′)=1(2π)3/2exp(−12r′Tr′)(29)The integration region V becomes V′ as (r′−L−1μ)TLTL(r′−L−1μ)≤r2(30)which is an ellipsoid. Then, μ′=L−1μ and A=LTL are the center and the quadratic matrix of the transformed ellipsoid, respectively. This process is illustrated in Fig. 1b.2. Rotation TransformA can be further diagonalized as QTAQ=Λ, where Q is an orthogonal matrix. Λ=diag{λ1,λ2,λ3}, where λ1, λ2, and λ3 are the three eigenvalues of A. Let r′=Qr′′ and then the PDF is subsequently transformed into p(r′′)=1(2π)3/2exp(−12r′′Tr′′)(31)The integration region now becomes (r′′−QTμ′)Tdiag{λ1r2,λ2r2,λ3r2}(r′′−QTμ′)≤1(32)which is still an ellipsoid with the following properties: The geometric center is μ′′=QTμ′.The lengths of the three semi-axes, designated as a′′, b′′, and c′′, respectively, are a′′=rλ1,b′′=rλ2,c′′=rλ3(33)The three symmetrical axes are parallel to the coordinate axes.This process is illustrated in Fig. 1c.3. Equivalent Volume ReplacementThe attention now turns to solving the integration of the standard 3-D Gaussian PDF over a general ellipsoid. As there is no explicit formula to solve this integral, three different approaches can be used to solve it: the accurate but time-consuming numerical integration; the fast but inaccurate analytical approximation; the semi-analytical method adopting numerical integration in some dimensions and the analytic approximation along the rest.As depicted in Fig. 1d, consider the externally tangent cuboid (the dashed red rectangle) of the ellipsoid described by Eq. (32). If we approximate the integration region as this cuboid, the approximated value will be larger than the actual value. This is because the volume of the externally tangent cuboid is larger than the volume of the ellipsoid. To mitigate this, we proportionally contract the three symmetrical axes of the cuboid and force its volume to be the same as the original ellipsoid. The contracted cuboid is designated as V′′′ (the solid blue rectangle).The volume of the ellipsoid is (4/3)πa′′b′′c′′, whereas the volume of the externally tangent cuboid is 8a′′b′′c′′. Denote a′′′, b′′′, and c′′′ as the half lengths of the three symmetrical axes of V′′′. Then, a′′′=π63a′′,b′′′=π63b′′,c′′′=π63c′′(34)The integration of Eq. (28) is now I≈∫∫∫V′′′p(r′′) dr′′=∫x0−a′′′x0+a′′′N(x;0,1) dx∫y0−b′′′y0+b′′′N(y;0,1) dy∫z0−c′′′z0+c′′′N(z;0,1) dz(35)where x0=μ′′(1), y0=μ′′(2), and z0=μ′′(3). It is noted that the result is decoupled along the three axes, which is one of the advantages of using the equivalent volume cuboid approximation.Let Ix, Iy, and Iz represent the integrations along the X, Y, and Z axes, respectively. Ix is Ix=∫x0−a′′′x0+a′′′N(x;0,1) dx(36)where N(x;0,1)=(1/2π)exp(−x2/2) is the standard univariate Gaussian PDF. The X-axis integration Ix can be explicitly expressed as Ix=∫x0−a′′′x0+a′′′N(x;0,1) dx=12[erf(x0+a′′′2)−erf(x0−a′′′2)](37)where erf(⋅) is the Gauss error function. The same results can be analogously obtained for Iy and Iz.Finally, the integral is approximated as I≈Ix⋅Iy⋅Iz(38)Substituting Eq. (38) into Eq. (25) for all the mixture components yields the final instantaneous collision probability at time t. Although the entire process is presented for a 3-D Gaussian PDF integral, it can also be directly applied to the 2-D case. Therefore, our method can be used to compute the cumulative collision probability with the linear collision method.V. Simulation ResultsA. Geosynchronous Equatorial Orbit CaseThe test case 4 presented in [1] is first used to evaluate the proposed method. This case involves two geosynchronous equatorial orbit (GEO) satellites with a small relative velocity (about 0.02 m/s). Because of the very small position and velocity errors at epoch, after about 3 days of propagation, both the primary and secondary target uncertainties should still remain Gaussian. Table 1 shows the initial and final position UCMs with 100,000 Monte Carlo particles. It is clear that the propagated uncertainty can be described well using a single Gaussian PDF.A one-on-one Monte Carlo simulation is used as a benchmark to verify the end-to-end fidelity of our approach. First, a sample set including nMC samples is randomly generated for each target according to its initial state PDF. These samples are then propagated simultaneously to the desired time via numerical integrations. Finally, the one-on-one distances between the two propagated sample sets are examined. If a pair of samples has a distance less than the radius of the hard body sphere, a collision is confirmed. The sample volume is determined according to [34] nMC>4(e−2)[(1−Pc)Pc]ε2Pc2log(21−γ)(39)where Pc is the collision probability, ε is the percentage error, and γ is the confidence level. A 1% error with 95% confidence level leads to a sample volume of 1E7 when Pc is assigned with 0.01.Figure 2 depicts the results for the instantaneous collision probability versus time with four different approaches. We find that the numerical integration results (denoted as "Numeric") agree well with the Monte Carlo results (denoted as "MC"). Compared with the result of Monte Carlo, the maximum collision probability by the equivalent cuboid method is overestimated approximately by 2%, whereas this probability by Chan's method is underestimated approximately by 86%. The proposed method (denoted as "Cuboid") has a better accuracy than Chan's approach (denoted as "Chan").Fig. 2 Instantaneous collision probability results for the GEO case.When applying the equivalent volume approximation in the isotropic probability space (as illustrated in Fig. 1d), both the volume loss and gain will occur. If the lost volume is more significant, the approximated Pc will be less than the exact value, and vice versa. In this case, Pc is underestimated around t=58 min. The reason is that one end of the transformed hard body ellipsoid is close to the origin. After applying the equivalent volume replacement, the corresponding face of the cuboid directly jumps over the origin, resulting in a dramatic loss of the collision probability.B. Low Earth Orbit CaseThe second test case consists of two low Earth orbit (LEO) satellites with orbit altitudes of about 400 km. For each satellite, the least-square orbit determination process is first conducted with simulated measurements from three radars to yield a more practical covariance matrix. At epoch, the position and velocity accuracies measured in 3-D root-mean-square error for the primary target are 0.965 m and 0.005 m/s, respectively, whereas those for the secondary target are 18.467 m and 0.088 m/s, respectively. The relative speed is approximately 4 km/s at the closest approach time, and the minimum separation is about 14 km. The hard body sphere radius is 200 m.The propagated uncertainty of