Nonlinear Probability Tool

The Nonlinear Probability Tool can be used to compute conjunction probability for nonlinear relative motion. The tool allows you to test for linear relative motion and provides two methods for computing conjunction probability: adjoining parallelepipeds (bundles) and adjoining tubes (cylinders).

If you have not already computed the close approach analysis for the selected Advanced CAT object, click Compute AdvCAT to perform that calculation before working with the tool. The Primary Object, Secondary Object, and Time of Closest Approach fields should be populated with valid selections for the selected Advanced CAT object. The fields beneath the object lists pertain to one or more of the functions that you can perform with the tool, as described in the sections below.

Testing for Linear Relative Motion

Since computing nonlinear probability requires very intensive calculations, it is usually a good idea to ensure that such a computation is necessary. Define a Fractional Limit of Probability (FLP) that is within accuracy requirements for your intended operations; FLP is the absolute difference between the linear and nonlinear results divided by the linear results. Click Compute Linearity Test to perform the test and review the results in the TCA Rel Val and Approx Min Rel Val fields to check whether a nonlinear computation is warranted.

The linearity test is comprised of two consecutive tests to find the minimum relative velocity needed at the time of closest approach (TCA) to ensure that a pre-specified probability difference will not be exceeded either forward or backward in time while the n- s shell is traversed. The first test is a coarse assessment that is performed only on the covariance growth using simple, two-body, orbital dynamics to approximate the bounds of linear motion based on the threshold FLP. This bounding is followed by a refined assessment using the full force models to predict relative orbital motion and covariance changes. The minimum relative velocity is then determined to ensure the Fractional Limit of Probability will not be exceeded. The linear motion assumption is valid if the encounter's relative velocity exceeds the minimum.

Coarse Assessment

For the bounding assessment, the relative motion is treated as linear but the covariance is allowed to grow. The assessment begins with the objects' positions, velocities, and covariance data at time of closest approach (TCA) where time is reset such that t=0. The examination region is bounded by specifying the size of the n- s shell that must be traversed in the End Sigma field. An initial estimate of time change Dt init is made by dividing the distance to traverse n- s by the magnitude of relative velocity. Beginning and end times (t start, t end) are preliminarily determined by using the relative velocity at time of closest approach and holding the n- s shell static to determine entry and exit. The next step is to compute the linear-motion probability at TCA, P(TCA). A plane perpendicular to the relative velocity vector is formed and the combined object and covariance ellipsoid are projected onto this encounter plane. The coarse approximation assumes linear relative motion coupled with a dynamic covariance such that the combined object center remains fixed at (x, y) and the 2x2 combined covariance C2 is the only parameter that changes with time.

Using the data at TCA, constant relative velocity is assumed and the relative position computed for t start. Each object's covariance is propagated to t start. This data is then used to find the associated probability at t start , P(t start). If the fractional tolerance between P(TCA) and P(t start) is exceeded then t start is assigned half its previous value and the process repeated until within the tolerance. When completed, t start becomes an initial estimate for how quickly the n- s shell must be traversed from TCA for the conjunction to be considered linear. The same process is repeated for t end.

Refined Assessment

Given P(t start), P(TCA), and P(t end) from the coarse assessment and realizing that the P curve is convex¹, a quadratic curve fit is employed to more closely associate the start and end times with the fractional probability tolerance. A new initial estimate for how quickly the n- s shell must be traversed from TCA is found by halving the difference between the fitted start and end times to produce Dt coarse. A coarse assessment of the minimum relative velocity needed to meet the fractional tolerance is determined by scaling the TCA relative velocity by Dt init/ Dt coarse. The secondary object's velocity is then modified to reflect this; the choice of modifying the primary or secondary is arbitrary. These new times and this new velocity are used with complex orbital force models to propagate positions and covariances from TCA to find the associated P(t start) and P(t end). A quadratic curve fit is again employed to refine the traversal time and produce Dt fine. The minimum relative velocity needed is recomputed by scaling the original relative velocity by Dt init/ Dt fine.

To do this, the unitized parameter t is computed from a quadratic curve fit as follows:

(1a)

(1b)

(2)

The coarse assessment data is inserted into the above equations to find Dt coarse:

(3)

The original relative velocity is scaled by Dt init/ Dt coarse. and the secondary object's velocity modified to reflect this. End and start times are reset to +/ - Dt coarse. P(t start) and P(t end) are recomputed from TCA using the secondary's modified velocity and complex orbital force models to propagate positions and covariances. Equations 1 and 2 are employed again to produce a new t which produces a refined Dt as follows:

(4)

The approximate, minimum, relative velocity needed to ensure enough linearity is finally determined by scaling the original relative velocity by Dt init/ Dt fine. This new velocity is consistent with the linear-motion assumptions based on the user-defined fractional probability threshold. If the encounter's relative velocity equals or exceeds this minimum, then the linear-motion assumption holds. If not, then a probability method that does not assume linear relative motion should be used. For more information please see Reference ².

Method of Adjoining Parallelepipeds (Bundles)

The collision tube is represented by bundles of abutting parallelepipeds. Each parallelepiped end is adjusted to form a compound miter where neighboring tubes meet, thereby reducing any gaps or overlaps. The approach applies to all relative motion and is coupled with a modified error-function method³ to allow any object shape. The method begins with object position, velocity, and covariance data at the time of closest approach. Propagation is done forward/backward in time until a limit you define is reached. For each time step the tube sections are sufficiently small enough so that, over the interval, the relative motion can be assumed linear and the covariance constant. The probability of each parallelepiped is computed and summed to get the probability of the tube section. All sections are summed to produce the encounter probability.

General Method

Geometric projections determine the end points of each parallelepiped. Let r1, r2, and r3 be three consecutive points along the relative trajectory in the VNC frame. Determine the unit vectors from r1 to r2 (axis12 for the first tube) and r2 to r3 (axis23 for the second tube). Rotate the axes to a new frame (denoted by suffix r) where the z component is aligned with axis12 such that after rotation axis12r is (0 0 1).

Define axis13r as the sum of axis12r and axis23r; the compound miter is perpendicular to axis13r and passes through r2r. In the new frame the r2r end point adjustment dz for each parallelepiped is found by examining the first tube's off-axis positions dx and dy through the equation:

(1)

The center of the parallelepiped's face is shifted from r2r by (dx dy dz), placing it on the surface of the compound miter.

At every time step the two-dimensional probability P 2d is computed by aligning the parallelepiped sides with the projected covariance axes. This eliminates covariance cross-correlation terms so that Equation 1 can be used for each of the two axes individually and the results multiplied to produce P 2d. The parallelepiped ends as adjusted by Equation 1 are transformed to the Mahalanobis space where Equation 2 is used to compute the long-axis probability P 1d:

(2)

Because of the approach taken to produce P 2d, this method can be extended to accommodate any complex object shape (convex, concave, spiral, hollow, etc.).

Sectional Computation

Three-dimensional position and velocity data of each object, as well as their 6x6 covariance matrices, are required with the assumption that all starting data are in the Earth Centered Inertial (ECI) frame then transformed to the primary object's Velocity-Normal-Co-Normal (VNC) frame where the relative distance vector is computed. Suitable incremental limits should be set for each time step with you specifying the computational stopping condition by time limit and/or encounter region. The computational algorithm is as follows:

Initially propagate all to Time of Closest Approach (TCA) in Earth Centered Inertial (ECI) frame

1. Convert starting data to VNC frame of primary object
2. Assign original relative position in VNC frame to r2
3. Determine relative position r1_ECI & covariance by propagating back one time step from TCA
4. Convert propagated data to VNC frame of primary object
5. Assign relative position in VNC frame to r1
6. Determine relative position r3_ECI & covariance by propagating forward one time step from TCA
7. Convert propagated data to VNC frame of primary object
8. Assign relative position in VNC frame to r3

Begin iteration

1. Propagate forward one time step from r3_ECI to determine relative position r4_ECI & covariance
2. Convert propagated data to VNC frame of primary object
3. Assign relative position in VNC frame to r4
4. Create unit vector from r1 to r2, label it axis12
5. Create unit vector from r2 to r3, label it axis23
6. Create unit vector from r3 to r4, label it axis34
7. Create vector from summation of axis12 and axis 23, label it axis13
8. Create vector from summation of axis23 and axis 34, label it axis24
9. Compute necessary rotation matrix to align new z component with relative velocity (axis23) while simultaneously decoupling new x and y components on projected covariance.
10. Rotate r2, r3, axis23, axis13, axis24, and 3x3 positional covariance (C3) associated with r2 to new frame while denoting rotated data with an r suffix (r2r, r3r, axis23r, axis13r, axis24r, C3r)
11. Compute necessary rotation/scaling matrix to go from new frame to Mahalanobis space where the z component is aligned with the relative velocity vector, label it T_maha
12. Middle tube axis endpoints are r2r and r3r: [xm, ym, zm2]=r2r & [xm, ym, zm3]=r3r
13. Find z component of tube's central axis ends using T_maha transformation, label them zm_start & zm_end
14. For each pixel of combined object determine its width, height, and off-axis central position (dx, dy)
15. Use r2r, axis13r, dx and dy to find dz2 to define one end of parallelepiped [xm+dx, ym+dy, zm2-dz2]
16. Find z component of parallelepiped end using T_maha transformation, label it z_start
17. Use r3r, axis24r, dx and dy to find dz3 to define other end of parallepiped [xm+dx, ym+dy, zm3-dz3]
18. Find z component of parallelepiped end using T_maha transformation, label it z_end
19. Find parallelepiped's 2D probability (face) centered at [xm+dx, ym+dy] using corresponding width and height
20. Find parallelepiped's 1D probability (length) using z_start and z_end
21. If sign(zm_end - zm_start) is opposite of sign(z_end - z_start) then there is overlap
22. Negate parallelepiped's 1D probability
23. Multiply 1D and 2D probabilities and add to running sum
24. Reassign r2 to r1, r3 to r2, r4 to r3; do likewise for covariances
25. Repeat until final limit reached (time, number of orbits, encounter shell limit, . . .)

This iterative procedure must be done twice, once forward in time from TCA and once backward in time.

Practical Considerations

Several intermediate limits must be addressed. The Max Time Step sets the sectional time step upper limit, the Max Angular Bend limits the angular difference between abutting cylinders, the Max Sigma Step limits individual cylinder length based on the traversed Mahalanobis distance. You must also specify how many parallelepipeds are needed to adequately represent the combined object space; the Resolution field defines this granularity for the two-dimensional probability computation by determining grid size. At the beginning of each iteration, the time step is increased by 25% to prevent small iterative steps required for one part of the trajectory from retarding calculations elsewhere. If the incremental results exceeded any limit, the time step is halved and the incremental calculations repeated. The Integration Duration field defines the span of time that will be evaluated - half of it forward from the TCA and half backward; by default, this span is defined as one half of the orbital period. Once the limit is reached, the entire process can be repeated with all incremental limits halved. If the results agreed to within two significant figures, the calculated probability should be considered suitable. The computation of parallelepipeds can be considerably longer than the adjoining tubes (cylinders) method. For more information please see Reference ⁴.

Method of Adjoining Tubes (Cylinders)

Linear relative motion methods for computing satellite collision probability are extended in this method to accommodate nonlinear relative motion in the presence of changing position and velocity uncertainties. For linear relative motion, the probability along the relative velocity vector (collision tube) is unity and is conveniently removed from the calculations. For nonlinear motion, that dimension is reintroduced. This is done by breaking the collision tube into small sections, computing the probability associated with each section, and then summing. In the absence of an external ephemeris file, the HPOP propagator is used to propagate satellite position, velocity, and associated covariances.

General Method

The general method of adjoining tubes begins with object position and velocity data at the time of closest approach. Propagation is done forward/backward in time until a limit you determine is reached. The limit can be based on a standard deviation threshold (such as 3 s, 6s, or 10 s) or a specified time (such as one quarter of an orbital period). For each time step the tube sections should be sufficiently small enough so that, over the interval, the relative motion can be assumed linear and the covariance constant. For each section, a two-dimensional probability is computed by projecting the combined object shape onto a plane perpendicular to the relative velocity. In addition, a one-dimensional probability is computed along the relative velocity vector by determining the component position from the mean at each end of the tube and then dividing by the standard deviation for that axis, thus producing each endpoint's Mahalanobis distance. The product of these probabilities yields the sectional probability. All sectional probabilities are summed until the time and/or sigma limit is reached. The probability of each cylinder is determined by multiplying the two-dimensional linear-motion probability by the sectional (relative velocity axis) probability.

The tubes have no gaps when dealing with linear relative motion. For such cases, the nonlinear results will perfectly match the linear-motion probability for constant covariance and objects of any shape. Nonlinear motion causes gaps and overlaps where the tube sections meet. If the relative motion track bends towards the covariance ellipsoid center, then the overlapping sections will occur in regions of greater probability density with the gaps occurring in regions of lesser probability density. Although the gap and overlapping volumes are almost equal, the resulting summation causes an over inflation of the probability. If the relative motion track bends away from the covariance ellipsoid center, then the probability for cylindrical tubes will be underestimated because the gap is in a region of higher probability density. The amount of error will vary based on the degree of bending/overlap relative to probability density.

There are several choices you should carefully consider when implementing this method. The limits and time step must be selected to ensure adequacy for the intended analysis. For a very large time limit and cyclical relative motion it is possible to retrace the same path through the covariance space. An example would be one satellite circling the other in formation for hundreds of revolutions. The collision tube would continually trace over itself; if care is not taken, the single revolution probability could be summed hundreds of times. To avoid this, it is suggested that the total time limit not exceed one-half of the primary vehicle's orbital period. A large intermediate step can also cause errors if the sectional motion is not sufficiently linear or the sectional covariance is not sufficiently constant. A simple test for sufficiency is to halve the intermediate step limits and repeat the analysis. If the probability differences are within your tolerance, then the intermediate steps are adequate.

Sectional Computation

Three-dimensional position and velocity data of each object, as well as their 6x6 covariance matrices, are required. At each time step all data is converted to the Velocity-Normal-Co-Normal (VNC) frame.

To compute the sectional probability of each tube, all data is propagated for the given time step. If the angular difference between the previous and current relative velocity vectors exceeds the angular limit, the time step is halved and this process repeated. A coordinate transformation is accomplished to align the z axis with the relative velocity vector. The one-dimensional, z-axis, Mahalanobis distances of the cylinder endpoints (M i, M f) are used to compute long-axis probability P 1d from

(1)

If the endpoint differences should exceed the maximum change in long-axis sigma then the time step is halved and this process repeated. Projection of the collision tube onto the encounter plane, defined by the x-y axes, produces the necessary information to generate two-dimensional collision probability. The one- and two-dimensional probabilities are then multiplied to determine the sectional collision probability. This process is repeated until either final limit is reached.

For this analysis, at every time step the objects, their positions, and their covariances are transformed to the VNC frame, where the relative distance vector is computed, then to the Mahalanobis space. For each tube section, a two-dimensional probability P 2d is computed by projecting the combined object shape onto a plane perpendicular to the relative velocity. The product of P 1d and P 2d yields the sectional probability of each tube. This requires accommodating the elliptical object footprint in the Mahalanobis encounter plane.

One of the challenges when working in Mahalanobis space is the combined object distortion that results from re-scaling the combined positional covariance to have unit variance along each dimension. The transformation is accomplished by determining the eigen values and associated eigenvectors of the combined positional covariance matrix. The eigen vectors are used to reorient all data while diagonalizing the covariance. All data is then scaled using the corresponding eigen values so that the covariance matrix becomes the identity matrix. The spherical object in Cartesian space becomes an ellipsoid, the circle an ellipse. The method of adjoining tubes requires a two-dimensional probability calculation for each segment, but in this space the object footprint may no longer be circular. Following a derivation like the one explained in Patera's 2005 Engineering Note⁵, a numerical series is derived to represent the combined object footprint as an ellipse in Mahalanobis space. For mathematical details refer to the Help section “Patera's Method (Modified).”

Practical Considerations

Several intermediate limits must be addressed. The Max Time Step sets the sectional time step upper limit, the Max Angular Bend limits the angular difference between abutting cylinders, the Max Sigma Step limits individual cylinder length based on the traversed Mahalanobis distance. At the beginning of each iteration, the time step is increased by 25% to prevent small iterative steps required for one part of the trajectory from retarding calculations elsewhere. If the incremental results exceeded any limit, the time step is halved and the incremental calculations repeated. The Integration Duration field defines the span of time that will be evaluated - half of it forward from the TCA and half backward; by default, this span is defined as one half of the orbital period. Once the limit is reached, the entire process can be repeated with all incremental limits halved. If the results agreed to within two significant figures, the calculated probability should be considered suitable. For more information please see Reference ².

Several intermediate limits must be addressed. The Max Time Step sets the sectional time step upper limit, the Max Angular Bend limits the angular difference between abutting cylinders, the Max Sigma Step limits individual cylinder length based on the traversed Mahalanobis distance. As an initial recommendation, set Max Time Step to 10 seconds, Max Angular Bend to one degree, and Max Sigma Step to 0.25 s. At the beginning of each iteration, the time step is increased by 25% to prevent small iterative steps required for one part of the trajectory from retarding calculations elsewhere. If the incremental results exceeded any limit, the time step is halved and the incremental calculations repeated. It is recommended that the time limit be initially set to one half of the orbital period (quarter orbit forward and quarter orbit backward from TCA) and the encounter shell's scale factor set to 10 (n=10). Once the limit is reached, the entire process can be repeated with all incremental limits halved. If the results agreed to within two significant figures, the calculated probability should be considered suitable. For more information please see Reference ⁶.