Method for calculating collision probability of rising section of aircraft in atmosphere based on direct method
Through the method based on the direct method, the probability of collision in the aircraft in the atmospheric atmosphere is calculated, and the problem of failure to consider the uncertainty of trajectory speed in the prior art is solved, and more accurate and efficient collision probability calculation is achieved, which is suitable for short-term and long-term encounters of aircraft.
Patent Information
- Application Number
- CN202510451232.7
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-04-11
- Publication Date
- 2025-07-11
AI Technical Summary
When calculating the probability of an aircraft's atmospheric collision, the prior art fails to effectively consider the velocity uncertainty of the trajectory, resulting in inaccurate calculation results and low calculation efficiency.
Using a direct method, by rotating the state mean and error covariance matrix under the aircraft ballistic system to the geocentric solid system, the mean and covariance matrix of the relative states of the two aircraft is calculated, and the Lejand polynomial is used to calculate the collision probability change rate under the spherical coordinate system, generating a collision probability change rate sequence, and finally normalization is performed to improve calculation accuracy and efficiency.
The calculation accuracy of collision probability in the short-term and long-term encounters of the aircraft is improved, the calculation efficiency is improved, and the real possibility of the aircraft collision and its time change information are reflected through the normalized collision probability.
Smart Images
Figure CN120296304A_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of aerospace engineering, and particularly to a method for calculating the collision probability of an in - atmosphere aircraft during the ascending stage based on the direct method. Background Technique
[0002] For the trajectory operation state of a ground - takeoff aircraft, due to the existence of internal deviations and external disturbances in the dynamic system, the actual flight state may deviate from the nominal state. During short - time, large - scale takeoff and flight processes, the distance between aircraft may be too close, resulting in mutual interference or even collision between aircraft, thus causing the flight mission to fail. Therefore, calculating the collision probability between aircraft can be used to calculate the takeoff time, flight path planning, etc., which is of great significance for the safe flight of aircraft.
[0003] The calculation of the collision probability is based on the relative motion state between two aircraft: a safety domain is established with the centroid of the primary star as the center, and a combined error ellipsoid of the relative position is established with the centroid of the secondary star as the center. At a certain moment, the probability that the relative position state of the aircraft falls within the safety domain is the collision probability. When analyzing these random errors by statistical methods, the error distributions of the aircraft position and velocity states are generally described by the covariance matrix or the error ellipsoid. At the same time, it is also necessary to define the safety area of the aircraft itself. Currently, most research uses a spherical shape, ignoring the influence of the aircraft attitude on the collision possibility.
[0004] Nowadays, the calculation and analysis of the collision probability of in - atmosphere aircraft are mostly based on the calculation and analysis on their encounter plane. For example, a fast algorithm for calculating the collision probability of space debris based on space compression and infinite series projects the relative position error of two aircraft into their encounter plane, simplifies the three - dimensional spherical volume integral to a two - dimensional circular domain integral, and approximates the integral as an infinite series form based on space compression and infinite series methods, giving an analytical calculation model for the collision probability. By performing numerical integration along the relative motion direction of the two aircraft to eliminate the one - dimension parallel to the velocity, the three - dimensional integral is simplified to two - dimensions. Considering the relative velocity error, the three - dimensional volume integral is transformed into the form of spherical surface integral and time integral, giving a method for calculating the collision probability of space targets under spacecraft dynamics.
[0005] Based on this, the present invention proposes a method for calculating the collision probability of an in - atmosphere aircraft during the ascending stage based on the direct method, which considers the velocity uncertainty of the trajectory, making the method proposed by the present invention applicable to the cases of short - term and long - term encounters of aircraft, so as to improve the accuracy of the calculation results. Summary of the Invention
[0006] The object of the present invention is to provide a method for calculating the collision probability of an in-atmosphere aircraft during the ascending stage based on the direct method. Through the definition of the direct method, in the integral calculation of the position error distribution of the flight trajectory, the velocity uncertainty of the trajectory is considered, making the method applicable to both short-term and long-term encounters of aircraft, and the results are more accurate; the triple integral is simplified to a double spherical surface integral, and through polynomial summation and approximation, a numerical integral solution is provided, which improves the calculation efficiency while ensuring the calculation accuracy; the normalized collision probability is used to simultaneously reflect the true value of the aircraft collision possibility and the information about its change over time.
[0007] To achieve the above object, the present invention provides a method for calculating the collision probability of an in-atmosphere aircraft during the ascending stage based on the direct method, including the following steps:
[0008] Step S1: Rotate the state mean and error covariance matrix in the aircraft ballistic system to the Earth-centered Earth-fixed (ECEF) system;
[0009] Step S2: Calculate the mean and covariance matrix of the relative state of the two aircraft in the ECEF coordinate system;
[0010] Step S3: Calculate the collision probability based on the direct method;
[0011] Step S4: Calculate the rate of change of the collision probability in the spherical coordinate system using Legendre polynomials;
[0012] Step S5: Generate a sequence of the rate of change of the collision probability based on the above process and calculate the cumulative collision probability sequence;
[0013] Step S6: Calculate the maximum value of the rate of change of the collision probability and the maximum value of the cumulative collision probability according to the sequence of the rate of change of the collision probability and the cumulative collision probability sequence, and normalize the cumulative collision probability to obtain the normalized collision probability.
[0014] Preferably, in step S1, when rotating the state mean and error covariance matrix in the aircraft ballistic system to the ECEF system, based on the three Euler angles ν E , σ E , θ E of the aircraft ballistic system relative to the ECEF system, the rotation matrix R d2e is obtained as follows:
[0015]
[0016] Preferably, the state mean m d in the aircraft ballistic system is as follows:
[0017]
[0018] Where r dis the mean value of the 3D position state in the ballistic coordinate system; v d is the mean value of the 3D velocity state in the ballistic coordinate system;
[0019] Then the mean value of the state of the aircraft in the Earth-centered Earth-fixed coordinate system m e , is as follows:
[0020]
[0021] v e = R d2e v d ;
[0022]
[0023] where, [x0, y0, z0] is the coordinate of the origin of the ballistic coordinate system in the Earth-centered Earth-fixed coordinate system; r e is the mean value of the 3D position state in the Earth-centered Earth-fixed coordinate system; v e is the mean value of the 3D velocity state in the Earth-centered Earth-fixed coordinate system.
[0024] Preferably, the covariance matrix P d of the state error of the aircraft in the ballistic coordinate system is as follows:
[0025]
[0026] where, P r_d represents the position state covariance matrix in the ballistic coordinate system; P v_d represents the velocity state covariance matrix in the ballistic coordinate system; P rv_d represents the covariance matrix related to position and velocity in the ballistic coordinate system; P vr_d represents the covariance matrix related to velocity and position in the ballistic coordinate system; P r_d , P v_d , P rv_d and P vr_d are all 3×3 matrices;
[0027] Then the error covariance matrices related to the position, velocity, and position-velocity of the aircraft in the Earth-centered Earth-fixed coordinate system are as follows:
[0028]
[0029] where, P r_e represents the position state covariance matrix in the Earth-centered Earth-fixed coordinate system; P v_e represents the velocity state covariance matrix in the Earth-centered Earth-fixed coordinate system; P rv_e represents the covariance matrix related to position and velocity in the Earth-centered Earth-fixed coordinate system; P vr_e represents the covariance matrix related to velocity and position in the Earth-centered Earth-fixed coordinate system.
[0030] Therefore, the error covariance matrix P in the Earth-Centered Earth-Fixed coordinate system e , is as follows:
[0031]
[0032] Preferably, in step S2, in the Earth-Centered Earth-Fixed coordinate system, calculate the mean and covariance matrix of the relative states of the two aircraft. The specific process is as follows:
[0033] Step S21: In the Earth-Centered Earth-Fixed coordinate system, based on the state means m a and m b of the two flight trajectories, calculate the relative state mean μ of the two aircraft a and b, as follows:
[0034] μ = m a - m b ;
[0035] Step S22: Based on the state error covariance matrices P a and P b of the two flight trajectories, calculate the relative error covariance matrix Σ, as follows:
[0036] Σ = P a + P b .
[0037] Preferably, in step S3, calculate the collision probability based on the direct method. The specific process is as follows:
[0038] Step S31: According to the definition of the collision probability change rate p c (t) of the direct method for the collision probability of two space objects A and B at time t, it is as follows:
[0039]
[0040] Among them, S represents the total surface area of the sphere; p(x t ) represents the probability density function of the relative state vector x t of object A and object B; v n represents the projection of the relative velocity on the spherical unit normal vector; v n ≤ 0 corresponds to the direction flowing into the spherical surface; v t represents the relative velocity vector of object A and object B at time t; represents the unit normal vector of the surface area differential element dS;
[0041] Step S32: Based on the probability density function of the Gaussian mixture model, the relative state distribution p(x t ) between object A and object B is as follows:
[0042] p(xt ) = p g (x t ; μ t , Σ t );
[0043] Among them, μ t and Σ t respectively represent the mean and covariance matrix of the relative states of the two aircraft at time t;
[0044] Step S33: Substitute the above formula into the collision probability change rate to obtain the collision probability change rate.
[0045] Preferably, μ t , Σ t are decomposed into the product of the probability density function of the relative position and the probability density function of the conditional relative velocity, as follows:
[0046]
[0047] Among them, μ r,t and μ v,t respectively represent the mean of the relative position state and the mean of the relative velocity state; Σ r,t and Σ v,t respectively represent the covariance matrix of the relative position state and the covariance matrix of the relative velocity state; Σ rv,t and Σ vr,t respectively represent the covariance matrix of the correlation between the relative position and the relative velocity, and the covariance matrix of the correlation between the relative velocity and the relative position;
[0048] Then the Gaussian probability density function form of the relative state distribution p(x t ) between object A and object B is decomposed into the product of Gaussian components, as follows:
[0049] p g (x t ; μ t , Σ t ) = p g (r t ; μ r,t , Σ r,t )p g (v t '; (μ v,t )', (Σ v,t )');
[0050] v t ' = v t - Σ vr,t (Σ r,t ) -1 r t ;
[0051] (μv,t )′ = μ v,t -Σ vr,t (Σ r,t ) -1 μ r,t ;
[0052] (Σ v,t )′ = Σ v,t -Σ vr,t (Σ r,t ) -1 Σ rv,t ;
[0053] where r t represents the relative distance vector between two objects at time t;
[0054] Then the relative state distribution p(x t ) between objects A and B is as follows:
[0055] p(x t ) = p g (r t ; μ r,t , Σ r,t )p g (v t ′; (μ v,t )′, (Σ v,t )′);
[0056] Therefore, the rate of change of the collision probability p c (t) is as follows:
[0057] p c (t) = ∫ S p g (r t ; μ r,t , Σ r,t )ν(r t )dS;
[0058]
[0059] Converting the spherical surface integral of the rate of change of the collision probability in the above formula to the spherical coordinate system, we have dS = R 2 cosθdθdφ, where R represents the radius of the safety domain sphere. Then the rate of change of the collision probability p c (t) is as follows:
[0060]
[0061] where erf(·) is the error function.
[0062] Preferably, in step S4, the change rate of the collision probability in the spherical coordinate system is calculated using Legendre polynomials, and points are taken and segmented in two angular directions to calculate the nodes φ i and θ j of the Legendre polynomials and their corresponding weight coefficients Therefore, the change rate p c (t) of the collision probability between the two flight trajectories at the current moment is as follows:
[0063]
[0064] where n ij represents the unit normal vector of the spherical microelement corresponding to the Legendre node and is as follows:
[0065]
[0066] Preferably, in step S5, based on the above process, a sequence of the change rates of the collision probability is generated, and a sequence of the cumulative collision probabilities P c (t n ) is calculated as follows:
[0067] P c (t0) = 0;
[0068]
[0069] where P c (t0) represents the cumulative collision probability at the initial time t0; the integration method for calculating the values of the sequence of the cumulative collision probabilities P c (t n ) requires that the time series does not need to be equally spaced.
[0070] Preferably, in step S6, based on the sequence of the change rates of the collision probability and the sequence of the cumulative collision probabilities, the maximum value of the change rate of the collision probability and the maximum value of the cumulative collision probability are calculated, and the cumulative collision probability is normalized to obtain the normalized collision probability as follows:
[0071]
[0072] Therefore, the present invention adopts the above-mentioned method for calculating the collision probability of an in-atmosphere aircraft during the ascending phase based on the direct method. Through the definition of the direct method, in the integral calculation of the position error distribution of the flight trajectory, the velocity uncertainty of the trajectory is considered, making the method applicable to the cases of short-term and long-term encounters of aircraft, and the results are more accurate. The triple integral is simplified to a double spherical surface integral, and a numerical integral solution is provided through polynomial summation and approximation, which improves the calculation efficiency while ensuring the calculation accuracy. The normalized collision probability is used to simultaneously reflect the true value of the aircraft collision possibility and the information about its change over time.
[0073] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Description of the Drawings
[0074] Figure 1 is a flowchart of a method for calculating the collision probability of an in-atmosphere aircraft during the ascending phase based on the direct method;
[0075] Figure 2 is the rate of change of the collision probability in the embodiment of the present invention;
[0076] Figure 3 is the trajectory of the two aircraft in the embodiment of the present invention;
[0077] Figure 4 is the cumulative collision probability in the embodiment of the present invention;
[0078] Figure 5 is the normalized collision probability in the embodiment of the present invention. Detailed Embodiment
[0079] The technical solution of the present invention will be further described below with reference to the accompanying drawings and embodiments.
[0080] As Figure 1 shown, a method for calculating the collision probability of an in-atmosphere aircraft during the ascending phase based on the direct method includes the following steps:
[0081] Step S1, rotate the state mean and the error covariance matrix in the aircraft ballistic coordinate system to the Earth-centered Earth-fixed coordinate system.
[0082] Step S11, based on the three Euler angles ν E , σ E , θ E of the aircraft ballistic coordinate system relative to the Earth-centered Earth-fixed coordinate system, obtain the rotation matrix R d2e , as follows:
[0083]
[0084] Step S12, the state mean m in the aircraft ballistic coordinate system d, as follows:
[0085]
[0086] Among them, r d is the 3D position state mean in the ballistic coordinate system; v d is the 3D velocity state mean in the ballistic coordinate system.
[0087] Then the state mean m of the aircraft in the Earth-centered Earth-fixed coordinate system e , as follows:
[0088]
[0089] v e = R d2e v d (4);
[0090]
[0091] Among them, [x0, y0, z0] is the coordinate of the origin of the ballistic coordinate system in the Earth-centered Earth-fixed coordinate system; r e is the 3D position state mean in the Earth-centered Earth-fixed coordinate system; v e is the 3D velocity state mean in the Earth-centered Earth-fixed coordinate system.
[0092] Step S13, the covariance matrix P of the state error of the aircraft in the ballistic coordinate system d , as follows:
[0093]
[0094] Among them, P r_d represents the position state covariance matrix in the ballistic coordinate system; P v_d represents the velocity state covariance matrix in the ballistic coordinate system; P rv_d represents the covariance matrix related to position and velocity in the ballistic coordinate system; P vr_d represents the covariance matrix related to velocity and position in the ballistic coordinate system; P r_d , P v_d , P rv_d and P vr_d are all 3×3 matrices.
[0095] Then the error covariance matrices of the position, velocity, and position-velocity correlation of the aircraft in the Earth-centered Earth-fixed coordinate system are as follows:
[0096]
[0097] Among them, P r_e represents the position state covariance matrix in the Earth-centered Earth-fixed coordinate system; P v_e represents the velocity state covariance matrix in the Earth-centered Earth-fixed coordinate system; Prv_e Denote the position and velocity correlation covariance matrix in the Earth-Centered Earth-Fixed (ECEF) coordinate system; P vr_e Denote the velocity and position correlation covariance matrix in the Earth-Centered Earth-Fixed (ECEF) coordinate system.
[0098] Therefore, the error covariance matrix P in the Earth-Centered Earth-Fixed coordinate system e is as follows:
[0099]
[0100] Step S2: Calculate the mean and covariance matrix of the relative states of the two aircraft in the Earth-Centered Earth-Fixed coordinate system.
[0101] Step S21: Based on the state means m a and m b of the two flight trajectories in the Earth-Centered Earth-Fixed coordinate system, calculate the relative state mean μ of the two aircraft a and b as follows:
[0102] μ = m a - m b (9);
[0103] Step S22: Based on the state error covariance matrices P a and P b of the two flight trajectories, calculate the relative error covariance matrix Σ as follows:
[0104] Σ = P a + P b (10).
[0105] Step S3: Calculate the collision probability based on the direct method.
[0106] In the definition of the direct method, assume that two space objects collide at time t, V is the volume of the safety domain sphere, and r t is the relative distance vector of the two objects at time t. If the radius of the safety domain sphere is R, when r t ∈V, that is, ||r t || ≤ R, then a collision occurs.
[0107] Step S31: According to the direct method of collision probability, the definition of the collision probability change rate p c (t) is as follows:
[0108]
[0109] where S represents the total surface area of the sphere; p(x t ) represents the probability density function of the relative state vector x t ; v n represents the projection of the relative velocity on the spherical unit normal vector; vn ≤0 corresponds to the direction flowing into the spherical surface; v t represents the relative velocity vector of two objects at time t; represents the unit normal vector of the surface area differential element dS.
[0110] Step S32, let x a (t) = x a,t 、x b (t) = x b,t respectively represent the state vectors of two spatial objects A and B, and x(t) = x t is the relative state vector between objects A and B, where x t = x b,t - x a,t . The probability density function p(x t ) of the relative state vector x t is as follows:
[0111] p(x t ) = ∫p(x a,t , x b,t )dx a,t (13);
[0112] where p(x a,t , x b,t ) represents the joint probability density function of the time-varying state vectors of object A and object B.
[0113] Assume that the state vectors of object A and object B are independent of each other, and their joint probability density function is the product of the probability density functions of individual state vectors, then there is:
[0114] p(x a,t , x b,t ) = p(x a,t )p(x b,t ) (14);
[0115] where p(x a,t ) represents the probability density function of the state vector of object A; p(x b,t ) represents the probability density function of the state vector of object B.
[0116] Based on the probability density function of the Gaussian mixture model, the probability density functions of the state vectors of object A and object B are given as follows:
[0117] p(x a,t ) = p g (x a,t ; m a,t , P a,t ) (15);
[0118] p(xb,t ) = p g (x b,t ; m b,t , P b,t ) (16);
[0119] Among them, the Gaussian probability density function is as follows:
[0120]
[0121] Among them, m and P respectively represent the mean and covariance matrix of the state error; n is the dimension of the uncertain parameter.
[0122] Therefore, the joint probability density function of the time-varying state vectors of object A and object B can be written as:
[0123] p(x a,t , x b,t ) = p(x a,t )p(x b,t ) = p g (x a,t ; m a,t , P a,t )p g (x b,t ; m b,t , P b,t )(18);
[0124] Step S33: Use x b,t = x a,t + x t to represent the state vector of object B, then there is:
[0125] ∫p g (x a,t ; m a,t , P a,t )p g (x a,t + x t ; m b,t , P b,t )dx a,t = p g (x t ; μ t , Σ t )(19);
[0126] Among them, μ t and Σ t respectively represent the mean and covariance matrix of the relative state of the two aircraft at time t.
[0127] Therefore, when the state distributions of object A and B are given by the Gaussian mixture probability density function, then the relative state distribution p(x between object A and B in formula (13)t ) can be written as:
[0128] p(x t ) = p g (x t ; μ t , Σ t ) (20);
[0129] Therefore, substituting the above formula into formula (11) to find the change rate of collision probability p c (t).
[0130] In the present invention, in order to simplify this integral calculation formula, μ t , Σ t is decomposed into the product of the probability density function of the relative position and the probability density function of the conditional relative velocity, as shown below:
[0131]
[0132] where μ r,t and μ v,t respectively represent the mean of the relative position state and the mean of the relative velocity state; Σ r,t and Σ v,t respectively represent the covariance matrix of the relative position state and the covariance matrix of the relative velocity state; Σ rv,t and Σ vr,t respectively represent the covariance matrix of the correlation between the relative position and the relative velocity, and the covariance matrix of the correlation between the relative velocity and the relative position.
[0133] Then the Gaussian probability density function form of the relative state distribution p(x t ) between objects A and B is decomposed into the product of Gaussian components, as shown below:
[0134] p g (x t ; μ t , Σ t ) = p g (r t ; μ r,t , Σ r,t )p g (v t ′; (μ v,t )′, (Σ v,t )′) (23);
[0135] v t ′ = v t - Σ vr,t (Σ r,t ) -1 r t (24);
[0136] (μ v,t )′ = μ v,t - Σ vr,t (Σ r,t ) -1 μ r,t (25);
[0137] (Σ v,t )′ = Σ v,t - Σ vr,t (Σ r,t ) -1 Σ rv,t (26);
[0138] where r t represents the relative distance vector between two objects at time t.
[0139] Substituting Equation (23) into Equation (20), the relative state distribution p(x t ) between objects A and B can be written in the following form:
[0140] p(x t ) = p g (r t ; μ r,t , Σ r,t )p g (v t ′; (μ v,t )′, (Σ v,t )′) (27);
[0141] Substituting Equation (27) into Equation (11), the rate of change of the collision probability p c (t) is as follows:
[0142] p c (t) = ∫ S p g (r t ; μ r,t , Σ r,t )ν(r t )dS (28);
[0143]
[0144] Converting the surface integral in Equation (28) to spherical coordinates, we have dS = R 2 cosθdθdφ, where R represents the radius of the safety domain sphere. Then the rate of change of the collision probability p c (t) is as follows:
[0145]
[0146]
[0147] Among them, erf(·) is the error function.
[0148] Step S4: Calculate the change rate of the collision probability in the spherical coordinate system by using Legendre polynomials.
[0149] Take points and divide them in two angular directions respectively, and calculate the nodes φ i and θ j of the Legendre polynomials and their corresponding weight coefficients Therefore, the change rate p c (t) of the collision probability between two flight trajectories at the current moment is as follows:
[0150]
[0151] Among them, n ij represents the unit normal vector of the spherical microelement corresponding to the Legendre node, as follows:
[0152]
[0153] Step S5: Repeat the above steps along the time series, calculate and generate a sequence of the change rate of the collision probability. According to the sequence of the change rate of the collision probability varying with time, calculate the cumulative collision probability sequence P c (t n ), as follows:
[0154] P c (t0) = 0 (38);
[0155]
[0156] Among them, P c (t0) represents the cumulative collision probability at the initial time t0.
[0157] In addition, the integration method for calculating the values of the cumulative collision probability sequence P c (t n ) requires that the time series does not have to be equally spaced.
[0158] Step S6: According to the sequence of the change rate of the collision probability and the cumulative collision probability sequence, calculate the maximum value of the change rate of the collision probability and the maximum value of the cumulative collision probability and normalize the cumulative collision probability to obtain the normalized collision probability
[0159]
[0160] Embodiment
[0161] Based on the above-mentioned method for calculating the collision probability during the ascending stage of an in-atmosphere aircraft based on the direct method proposed by the present invention, this embodiment further illustrates the invention content in conjunction with the accompanying drawings.
[0162] First, for the initial boost phase trajectory (duration 60 s) of launching an aircraft from the ground, a collision probability prediction is carried out. The initial masses of both aircraft are set to m0 = 2500 kg, and the initial states of the trajectory are shown in Table 1.
[0163] Table 1 Initial states of the aircraft
[0164]
[0165] Set the prediction step size to 0.01 s, and the initial state errors of the two aircraft are as follows:
[0166] Position error standard deviation 0.1 m, speed magnitude error standard deviation 0.01 m / s, launch inclination error standard deviation 2°, launch declination error standard deviation 2°.
[0167] From the above initial values, the mean and covariance matrix of the error distribution of the trajectory state in the aircraft ballistic system are obtained. Based on the three Euler angles of the aircraft ballistic system relative to the Earth-centered Earth-fixed system, a rotation matrix is obtained, and the state mean and error covariance matrix in the aircraft ballistic system are rotated to the Earth-centered Earth-fixed system. In the Earth-centered Earth-fixed coordinate system, the mean and covariance matrix of the relative state of the two aircraft are calculated.
[0168] Based on the direct method to calculate the collision probability, set the radius R of the safety domain sphere to 30 m. μ t and Σ t are respectively the mean and covariance matrix of the relative state of the two aircraft at time t. μ t , Σ t are decomposed into the product of the probability density function of the relative position and the probability density function of the conditional relative velocity to obtain the collision probability change rate at time t.
[0169] Using Legendre polynomials to approximately calculate the spherical surface integral, point sampling and segmentation are carried out in two angular directions respectively, the nodes of the Legendre polynomials and their corresponding weight coefficients are calculated, and the collision probability change rate between the two flight trajectories at the current moment is obtained.
[0170] Repeat the above process along the time series, calculate and generate a sequence of collision probability change rates. According to the sequence of the collision probability change rate varying with time, calculate the cumulative collision probability sequence. The collision probability change rate reflects the trend of the collision probability changing with time, and the cumulative collision probability represents the real collision possibility. Therefore, through normalization, the information contained in the two parameters is reflected simultaneously. Based on this, the nominal launch trajectory and the collision probability results are obtained, as Figures 2 - 5 shown.
[0171] As Figure 2 shown, the maximum change rate of the collision probability is 7.77%, and the corresponding time is 50.65 s, which is the moment when the collision possibility is the greatest.
[0172] As Figure 3 shown, the intersection of the trajectories of the two aircraft is approximately 50.68 s. Therefore, the moment when the collision probability is the greatest is the moment of trajectory intersection.
[0173] As Figure 4 shown, the maximum cumulative collision probability is 12.73%, and this value is the true value of the collision probability. Comparing it with the Monte Carlo method, when 5000 sample points are selected, the maximum cumulative collision probability result based on the Monte Carlo method is 12.84%, and the absolute error between the two is 0.11%.
[0174] It should be noted here that in the collision warning of spacecraft orbits, according to the internationally commonly used collision probability threshold, the yellow warning threshold is 10 -5 , and the red warning threshold is 10 -4 . That is, when the collision probability is higher than 10 -4 , it means that a collision is very likely to occur.
[0175] Therefore, as Figure 5 shown, the above collision probability prediction results show that the collision possibility of the two aircraft is very high, and it is necessary to change the nominal trajectory to achieve avoidance, which is consistent with the actual results.
[0176] Therefore, the present invention adopts the above-mentioned method for calculating the collision probability of the ascending section of an aircraft in the atmosphere based on the direct method. Through the definition of the direct method, in the integral calculation of the position error distribution of the flight trajectory, the velocity uncertainty of the trajectory is considered, making the method applicable to the short-term and long-term encounter situations of the aircraft, and the results are more accurate; the triple integral is simplified to a double spherical surface integral, and through polynomial and approximation, a numerical integral solution is provided, which improves the calculation efficiency while ensuring the calculation accuracy; the normalized collision probability is used to simultaneously reflect the true value of the collision possibility of the aircraft and the information of its change over time.
[0177] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions of the present invention or make equivalent replacements, and these modifications or equivalent replacements do not make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.
Claims
1. A method for calculating the collision probability of an in - atmosphere aircraft during the ascent phase based on the direct method, characterized in that, Including the following steps: Step S1: Rotate the state mean and error covariance matrix in the aircraft ballistic system to the Earth-centered Earth-fixed (ECEF) system; Step S2: Calculate the mean and covariance matrix of the relative state of the two aircraft in the ECEF coordinate system; Step S3: Calculate the collision probability based on the direct method; Step S4: Calculate the rate of change of the collision probability in the spherical coordinate system using Legendre polynomials; Step S5: Generate a sequence of the rate of change of the collision probability based on the above process and calculate the cumulative collision probability sequence; Step S6: Calculate the maximum value of the rate of change of the collision probability and the maximum value of the cumulative collision probability according to the sequence of the rate of change of the collision probability and the cumulative collision probability sequence, and normalize the cumulative collision probability to obtain the normalized collision probability.
2. The method for calculating the collision probability during the ascending stage of an in-atmosphere aircraft based on the direct method according to claim 1, wherein In step S1, the state mean and error covariance matrix in the vehicle ballistic coordinate system are rotated to the Earth-centered Earth-fixed coordinate system. Based on the three Euler angles ν E , σ E , θ E of the vehicle ballistic coordinate system relative to the Earth-centered Earth-fixed coordinate system, the rotation matrix R d2e is obtained as follows:
3. A method for calculating the collision probability of an in-atmosphere aircraft during the ascending stage based on the direct method according to claim 2, characterized in that, The state mean value m in the flight vehicle ballistic system d , as shown below: where r d is the mean value of the 3D position state in the ballistic system; v d is the mean value of the 3D velocity state in the ballistic system; Then the state mean m of the aircraft in the Earth-centered Earth-fixed coordinate system e , is as follows: v e = R d2e v d ; where [x0, y0, z0] are the coordinates of the origin of the ballistic coordinate system in the Earth-centered Earth-fixed coordinate system; r e is the mean value of the 3D position state in the Earth-centered Earth-fixed system; v e is the mean value of the 3D velocity state in the Earth-centered Earth-fixed system.
4. A method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 2, wherein Covariance matrix \(P\) of the state error in the flight vehicle ballistic system d , as follows: Among them, P r_d represents the position state covariance matrix in the ballistic system; P v_d represents the velocity state covariance matrix in the ballistic system; P rv_d represents the covariance matrix of position and velocity in the ballistic system; P vr_d represents the covariance matrix of velocity and position in the ballistic system; P r_d , P v_d , P rv_d and P vr_d are all 3×3 matrices; Then the position, velocity, and error covariance matrix related to position and velocity of the aircraft in the ECEF coordinate system are as follows: Among them, P r_e represents the position state covariance matrix in the Earth-Centered Earth-Fixed (ECEF) coordinate system; P v_e represents the velocity state covariance matrix in the ECEF coordinate system; P rv_e represents the covariance matrix related to position and velocity in the ECEF coordinate system; P vr_e represents the covariance matrix related to velocity and position in the ECEF coordinate system; Therefore, the error covariance matrix P in the Earth-centered Earth-fixed coordinate system e is as follows:
5. A method for calculating the collision probability of an atmospheric aircraft during the ascending stage based on the direct method according to claim 1, characterized in that, In Step S2, in the ECEF coordinate system, calculate the mean and covariance matrix of the relative state of the two aircraft. The specific process is as follows: Step S21: In the Earth-centered Earth-fixed coordinate system, based on the state means m a and m b of two flight trajectories, calculate the relative state mean μ of two aircraft a and b as follows: μ = m a -m b ; Step S22: Based on the state error covariance matrices P a and P b , calculate the relative error covariance matrix Σ as follows: Σ = P a + P b 。 6. The method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 1, wherein In Step S3, calculate the collision probability based on the direct method. The specific process is as follows: Step S31. According to the direct method of collision probability, the rate of change p c (t) of the collision probability between two space objects A and B at time t is defined as follows: Among them, S represents the total surface area of the sphere; p(x t ) represents the probability density function of the relative state vector x of object A and object B t ; v n represents the projection of the relative velocity on the unit normal vector of the spherical surface; v n ≤ 0 corresponds to the direction flowing into the spherical surface; v t represents the relative velocity vector of object A and object B at time t; represents the unit normal vector of the surface area differential element dS; Step S32. Based on the probability density function of the Gaussian mixture model, the relative state distribution p(x t ) between object A and object B is as follows: p(x t ) = p g (x t ; μ t , Σ t ); where μ t and Σ t represent the mean and covariance matrix of the relative states of the two aircraft at time t, respectively; Step S33: Substitute the above formula into the rate of change of the collision probability to obtain the rate of change of the collision probability.
7. A method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 6, characterized in that, Decompose μ t , Σ t into the product of the probability density function of the relative position and the probability density function of the conditional relative velocity as follows: where, μ r,t and μ v,t represent the mean of the relative position state and the mean of the relative velocity state respectively; Σ r,t and Σ v,t represent the covariance matrix of the relative position state and the covariance matrix of the relative velocity state respectively; Σ rv,t and Σ vr,t represent the covariance matrix of the correlation between the relative position and the relative velocity, and the covariance matrix of the correlation between the relative velocity and the relative position respectively; Then the Gaussian probability density function of the relative state distribution p(x t ) between object A and object B is factorized into a product of Gaussian components as follows: p g (x t ; μ t , Σ t ) = p g (r t ; μ r,t , Σ r,t )p g (v t '; (μ v,t )', (Σ v,t )) v t ′ = v t -Σ vr,t (Σ r,t ) -1 r t ; (μ v,t )′ = μ v,t - Σ vr,t (Σ r,t ) -1 μ r,t ; (Σ v,t )′ = Σ v,t - Σ vr,t (Σ r,t ) -1 Σ rv,t ; where r t represents the relative distance vector of two objects at time t; Then the relative state distribution p(x t ) between objects A and B is as follows: p(x t ) = p g (r t ; μ r,t , Σ r,t ) p g (v t ′; (μ v,t )′, (Σ v,t )′); Therefore, the rate of change of the collision probability p c (t) is as follows: p c (t) = ∫ S p g (r t ; μ r,t , Σ r,t ) ν(r t ) dS; Convert the spherical surface integral of the collision probability change rate in the above formula to the spherical coordinate system, and we have dS = R 2 cosθdθdφ, where R represents the radius of the safety domain sphere, then the collision probability change rate p c (t) is as follows: where erf(·) is the error function.
8. A method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 1, characterized in that, In step S4, the change rate of the collision probability in the spherical coordinate system is calculated using Legendre polynomials. Sampling and segmentation are respectively carried out in two angular directions to calculate the nodes φ i , θ j of the Legendre polynomials and their corresponding weight coefficients Therefore, the change rate p c (t) of the collision probability between the two flight trajectories at the current moment is as follows: where n ij represents the unit normal vector of the spherical differential element corresponding to the Legendre nodes, as follows:
9. A method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 1, wherein, In step S5, based on the above process, a sequence of collision probability change rates is generated, and a cumulative collision probability sequence P c (t n ) is calculated as follows: P c (t0) = 0; Among them, P c (t0) represents the cumulative collision probability at the initial time t0; calculate the cumulative collision probability sequence P c (t n ) The numerical integration method requires that the time series does not need to be equally spaced.
10. A method for calculating the collision probability of the ascending section of an in-atmosphere aircraft based on the direct method according to claim 1, characterized in that, In step S6, according to the sequence of collision probability change rates and the sequence of cumulative collision probabilities, calculate the maximum value of the collision probability change rate and the maximum value of the cumulative collision probability and normalize the cumulative collision probability to obtain the normalized collision probability as follows: