CN113343436A - Method and system for calculating collision probability of Gaussian mixture covariance evolution - Google Patents

Method and system for calculating collision probability of Gaussian mixture covariance evolution Download PDF

Info

Publication number
CN113343436A
CN113343436A CN202110550985.5A CN202110550985A CN113343436A CN 113343436 A CN113343436 A CN 113343436A CN 202110550985 A CN202110550985 A CN 202110550985A CN 113343436 A CN113343436 A CN 113343436A
Authority
CN
China
Prior art keywords
gaussian
covariance
evolution
gaussian mixture
orbit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110550985.5A
Other languages
Chinese (zh)
Other versions
CN113343436B (en
Inventor
闫瑞东
龚建村
罗冰显
刘四清
王荣兰
师立勤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
National Space Science Center of CAS
Original Assignee
National Space Science Center of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by National Space Science Center of CAS filed Critical National Space Science Center of CAS
Priority to CN202110550985.5A priority Critical patent/CN113343436B/en
Publication of CN113343436A publication Critical patent/CN113343436A/en
Application granted granted Critical
Publication of CN113343436B publication Critical patent/CN113343436B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention belongs to the technical field of spaceflight, and relates to a collision probability calculation method and a collision probability calculation system for Gaussian mixture covariance evolution, wherein the method comprises the following steps: performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit; carrying out covariance evolution on each Gaussian unit to obtain corresponding orbital evolution error distribution represented by Gaussian mixture; and calculating collision probability among Gaussian units aiming at the orbit evolution error distribution represented by the Gaussian mixture, and performing weighted summation to obtain the total collision probability. Compared with the prior art, the method provided by the invention can improve the calculation precision of the collision probability no matter whether the orbit uncertainty presents Gaussian distribution or non-Gaussian distribution.

Description

Method and system for calculating collision probability of Gaussian mixture covariance evolution
Technical Field
The invention belongs to the technical field of spaceflight, and relates to a collision probability calculation method and system for Gaussian mixture covariance evolution.
Background
In the collision early warning engineering, the time when two space objects are closest to each other in the next several days needs to be predicted, the track state and the covariance need to be forecasted, and the collision probability needs to be calculated. Various collision probability calculation methods assume that the orbit error follows Gaussian distribution, so a linear method is generally adopted to forecast the orbit covariance, the method adopts first-order linear approximation to a nonlinear orbit dynamic system, and the influence of high-order terms is ignored. In fact, when the initial orbit uncertainty is large or the orbit prediction period is long, the predicted orbit covariance will no longer obey the gaussian distribution. Therefore, a method for calculating the orbit covariance evolution and the collision probability aiming at a large initial orbit covariance or a long orbit prediction period is needed to improve the calculation accuracy of the collision probability and reduce the false alarm rate of collision early warning.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides a collision probability calculation method for Gaussian mixture covariance evolution.
In order to achieve the above object, the present invention provides a collision probability calculation method of gaussian mixture covariance evolution, including:
performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit;
carrying out covariance evolution on each Gaussian unit to obtain corresponding orbital evolution error distribution represented by Gaussian mixture;
and calculating collision probability among Gaussian units aiming at the orbit evolution error distribution represented by the Gaussian mixture, and performing weighted summation to obtain the total collision probability.
As an improvement of the above method, the two spatial objects to be predicted are subjected to initial orbit covariance fitting based on a gaussian mixture algorithm, respectively, to obtain a gaussian mixture probability density function of each spatial object with respect to a gaussian unit; the method specifically comprises the following steps:
respectively calculating Gaussian probability density of each space object relative to the initial orbit based on the divided Gaussian units; fitting the covariance based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function P (x) of the space object, wherein the Gaussian mixture probability density function P (x) is as follows:
Figure BDA0003075387070000021
where N is the total number of Gaussian units, pgi(x；mi,Pi) Is the probability representation of the ith Gaussian unit, x represents the state quantity of the initial orbit and comprises a three-dimensional velocity vector and a three-dimensional position vector, miIs the mean value of the probability density of the ith Gaussian cell, PiIs the covariance of the probability density of the ith Gaussian cell, alphaiIs the weighted value of the probability density function of the ith Gaussian cell.
As an improvement of the above method, the performing covariance evolution on each gaussian unit obtains an orbit evolution error distribution represented by a corresponding gaussian mixture; the method specifically comprises the following steps:
and carrying out covariance evolution on each Gaussian unit of each space object by adopting an Unscented transformation method, taking the estimated value and covariance of the state quantity x of the initial orbit at the 0 th moment as initial values, taking a discretization nonlinear equation as a model, and carrying out Unscented transformation until the estimated value and covariance at the T th moment are obtained, so that propagation of the mean value and covariance of the Gaussian units is realized to obtain the orbit evolution error distribution expressed by Gaussian mixture.
As an improvement of the above method, the estimated value and covariance of the state quantity x of the initial orbit at the 0 th time are used as initial values, a discretization nonlinear equation is used as a model, and the Unscented transformation is performed until the estimated value and covariance at the T th time are obtained; the method specifically comprises the following steps:
according to the estimated value of the state quantity x of the ith Gaussian unit initial orbit at the k-1 moment
Figure BDA0003075387070000022
Sum covariance
Figure BDA0003075387070000023
1≤k≤T
Figure BDA0003075387070000024
Figure BDA0003075387070000025
Build Sigma dots, get the following formula:
Figure BDA0003075387070000026
Figure BDA0003075387070000027
Figure BDA0003075387070000028
wherein σk-1Is and
Figure BDA0003075387070000029
a column vector of the same type is used,
Figure BDA00030753870700000210
the square root of the matrix is represented, lambda is an adjustable scale parameter and is used for improving approximation accuracy, and n is a state quantity dimension; sigma point χ is obtained by Cholesky decomposition or eigenvalue decomposition of covariance matrixk-1(j) Wherein Cholesky decomposition requires that the covariance matrix is positive, j represents the jth column of the covariance matrix;
directly carrying out state and observation prediction according to the calculated Sigma point according to the discretization nonlinear model to obtain a predicted symmetric distribution Sigma sample point chik(j) (j ═ 0, … 2n) is:
χk(j)＝f[χk-1(j),k]
wherein f [. cndot. ] represents a time-dependent nonlinear orbit prediction function, and k represents the kth moment;
according to Sigma sample point χk(j) Calculating mean of states
Figure BDA0003075387070000031
Sum variance
Figure BDA0003075387070000032
Figure BDA0003075387070000033
Figure BDA0003075387070000034
Wherein, χk(0) Represents the 0 th column of the Sigma sample point χ k (j), alpha and beta are constant values, and the value range of alpha is more than or equal to 1 and less than or equal to 10-4In the case of the gaussian distribution, β is 2.
As an improvement of the above method, the trajectory evolution error distribution expressed by the gaussian mixture calculates collision probabilities through gaussian units, and performs weighted summation to obtain a total collision probability; the method specifically comprises the following steps:
the time interval t is obtained according to the following formula0→tfProbability of collision of fcComprises the following steps:
fc＝P0+P1
wherein, P0Is t0Initial probability of collision at time, P1Is t0＜t≤tfA set of collision probabilities for the time intervals;
P0expressed as the integral of the velocity along the sphere coordinate system, satisfies the following equation:
Figure BDA0003075387070000035
wherein p isg(r0；μr(t0),Pr(t0) ) mean value of μr(t0) Covariance of Pr(t0) A multi-dimensional probability density distribution function, subscript r denotes the spherical radius, r0The radius of the spherical integral is represented, R represents the combined radius of the two intersecting objects, and theta and phi respectively represent the azimuth angle and the elevation angle of the spherical integral. When the distance between the space objects P and S exceeds a threshold value0＝0；
P1Satisfies the following formula:
Figure BDA0003075387070000041
where θ, φ represents the azimuth and elevation angles of the spherical integral, respectively, v (r ', t) represents the relative integral velocity, r' represents the integral radius, satisfying the following equation:
Figure BDA0003075387070000042
wherein the content of the first and second substances,
Figure BDA0003075387070000047
is an expression of the original variable, and satisfies the following formula:
Figure BDA0003075387070000043
Figure BDA0003075387070000044
σ2＝r’T(Pv-PvrPr -1Pvr T)r’
wherein, PvCovariance P representing the ellipsoid of the joint errorrelThe 4 th to 6 th column parts, P, corresponding to the 4 th to 6 th rowsvrCovariance P representing the ellipsoid of the joint errorrelColumn 1 to column 3 portions corresponding to row 4 to row 6; η represents an integral variable;
a Lebedev integration method is adopted for spherical integration, and the integration changing along with time adopts a standard integration formula;
the total collision probability P is obtained according to the following formulacComprises the following steps:
Figure BDA0003075387070000045
wherein f isc(. is the probability of collision between two Gaussian mixture units, NP,NSRespectively representing the number of Gaussian units of the space object P and the space object S, i and j respectively representing the ith Gaussian unit of the space object P and the jth Gaussian unit of the space object S,
Figure BDA0003075387070000046
respectively representing the weighting coefficient of the ith Gaussian unit of the spatial object P, the weighting coefficient of the jth Gaussian unit of the spatial object S, and TCA (i, j) representing the intersection time of the ith Gaussian unit and the jth Gaussian unit.
A collision probability calculation system of gaussian mixture covariance evolution, the system comprising: the system comprises a Gaussian mixture method orbit covariance fitting module, a covariance evolution module and a collision probability calculation module; wherein the content of the first and second substances,
the Gaussian mixture method orbit covariance fitting module is used for performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit;
the covariance evolution module is used for carrying out covariance evolution on each Gaussian unit to obtain corresponding orbit evolution error distribution expressed by Gaussian mixture;
and the collision probability calculation module is used for calculating collision probability among the Gaussian units according to the orbit evolution error distribution expressed by the Gaussian mixture, and performing weighted summation to obtain the total collision probability.
A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the method when executing the computer program.
A computer-readable storage medium, in which a computer program is stored which, when executed by a processor, causes the processor to carry out the above-mentioned method.
Compared with the prior art, the invention has the advantages that:
compared with the prior art, the method adopts the Gaussian mixture algorithm to calculate the collision probability to reach the collision probability precision of Monte Carlo calculation no matter whether the uncertainty of the orbit presents Gaussian distribution or non-Gaussian distribution, and the Gaussian mixture covariance calculation method improves the collision probability calculation precision.
Drawings
FIG. 1 is a flow chart of a collision probability calculation method of Gaussian mixture covariance evolution according to an embodiment of the present invention;
FIG. 2 is a 3 σ error contour plot of a radial along-the-track planar projection of Monte Carlo samples and linear, unscented Kalman Filter, and Gaussian mixture methods.
Detailed Description
The method is realized by the following technical scheme:
a collision probability calculation method of Gaussian mixture covariance evolution is basically implemented as follows:
firstly, fitting the initial orbit covariance of the orbit by adopting a Gaussian mixture algorithm, and decomposing the initial orbit covariance into smaller Gaussian units; then carrying out covariance evolution on each Gaussian mixture unit by adopting an Unscented transformation method to obtain orbital evolution error distribution represented by Gaussian mixture; and finally, calculating collision probability among Gaussian mixture units aiming at the orbit error represented by the Gaussian mixture, and performing weighted summation to calculate the total collision probability.
The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and examples.
Example 1
As shown in fig. 1, embodiment 1 of the present invention provides a collision probability calculation method of gaussian mixture covariance evolution, which includes the following specific steps:
the method comprises the following steps that firstly, a nonlinear probability density function is approximated through weighted summation of a limited number of Gaussian probability density functions, and weighted values of different Gaussian units are determined by a numerical optimization method. The gaussian mixture probability density function is expressed as:
Figure BDA0003075387070000061
where N is the total number of Gaussian probability density unit functions, αiIs a weighted value, satisfies:
Figure BDA0003075387070000062
the mean and covariance of each term in the gaussian mixture function are evolved to quantify the error of the nonlinear system. Gaussian mixing is used in conjunction with other error propagation methods, and evolution of the mean and covariance is achieved by UT transform methods.
The method specifically comprises the following steps:
when the covariance of Gaussian mixture evolves, firstly a square root factor S is found to make SSTP. Then divide the square root factor into columns, let skIs the kth column of S. The execution univariate decomposition library is assigned a weight value
Figure BDA0003075387070000063
Mean value
Figure BDA0003075387070000064
And standard deviation of
Figure BDA0003075387070000065
The value of (c). When decomposing along the k-th axis of the square root factor, the component weights, mean, and covariance are:
Figure BDA0003075387070000066
Figure BDA0003075387070000067
Figure BDA0003075387070000068
wherein
Figure BDA0003075387070000069
By selecting a square root factor composed of a particular eigenvalue and eigenvector. Thus, the decomposition of the covariance matrix P is:
P＝VΛVT
the square root factor is determined as S ═ V Λ1/2Where λ is the diagonal matrix. Using matrix decomposition, a one-dimensional Gaussian mixture library edge one is selected for useThe feature vectors are decomposed. When decomposition is performed along the k-th axis of matrix decomposition, the component weights, means, and covariances used become:
Figure BDA0003075387070000071
Figure BDA0003075387070000072
Pi＝VΛiVT
wherein v iskIs the kth feature vector of P. LambdaiIs the ith set of eigenvalues for the new component:
Figure BDA00030753870700000715
step two, adopting the Unscented transformation to estimate the system state x at the moment k-1
Figure BDA0003075387070000073
Sum covariance
Figure BDA0003075387070000074
And taking a discretization nonlinear equation as a model, and performing Unscented transformation to realize the propagation of the mean value and the covariance of the Gaussian mixture unit. Wherein
Figure BDA0003075387070000075
And
Figure BDA0003075387070000076
mean and covariance of the gaussian mixture unit.
Figure BDA0003075387070000077
Figure BDA0003075387070000078
The method specifically comprises the following steps:
(1) construction of Sigma dots
Figure BDA0003075387070000079
Figure BDA00030753870700000710
Figure BDA00030753870700000711
In the formula (I), the compound is shown in the specification,
Figure BDA00030753870700000712
for the square root of the matrix, the convention is
Figure BDA00030753870700000713
The covariance matrix is determined by Cholesky decomposition or eigenvalue decomposition (singular value decomposition) of the covariance matrix, wherein Cholesky decomposition requires that the covariance matrix be positive. The index j indicates the jth column of the matrix and is
Figure BDA00030753870700000714
The column vector of the same type takes values from 0 to 2 n. And lambda is an adjustable scale parameter, the approximation precision can be improved by adjusting the parameter, and n is the number of estimators.
(2) Nonlinear propagation with Sigma spots
And directly carrying out state and observation prediction according to the calculated Sigma point according to a discretization nonlinear model, wherein the correspondingly generated Sigma sample points are as follows:
χk(j)＝f[χk-1(j),k]
middle X typek(j) (j ═ 0, … 2n) is the predicted Sigma point of the symmetric distribution.
(3) Calculating state mean and variance from sample points
Figure BDA0003075387070000081
Figure BDA0003075387070000082
Step three, calculating the collision probability by a Gaussian mixture method, namely calculating the collision probability between two target weighted gaussians
The method specifically comprises the following steps:
the basic formula for calculating the collision probability by a Gaussian mixture method is as follows:
Figure BDA0003075387070000083
fc(. cndot.) is the probability of collision between two Gaussian mixture units. When the two objects are closest, the relative speed of the two objects is perpendicular to the relative position of the two objects. The two objects meet in a plane perpendicular to the relative velocity, which is called the meeting plane. And (4) projecting the covariance information of the two targets into the intersection plane, and converting the three-dimensional collision probability calculation problem into two-dimensional collision probability calculation. And establishing the joint error distribution and the joint relative size of the two targets in the intersection plane, and calculating the collision probability through numerical integration. The collision probability of each target's mixed gaussian covariance relative to the covariance of the other target is calculated in the above equation,
Figure BDA0003075387070000085
and
Figure BDA0003075387070000086
respectively, the weighting coefficients of the neutron gaussians in the two targets.
Calculating the change rate of the instantaneous collision probability by integrating on the combined error ellipsoid of the two targets, and calculating the collision probability f by carrying out numerical integration in the intersection time interval of the two targetsc(. cndot.). The central point of the combined error ellipsoid is located at the second target position, and the relative state quantity x thereofrelRelative position rrelRelative velocity vrelCovariance of the combined error ellipsoid is Prel。
xrel＝xp-xs
rrel＝rp-rs
vrel＝vp-vs
Prel＝Pp-Ps
Figure BDA0003075387070000084
At a time interval t0→tfThe collision probability of (1) is:
fc＝P0+P1
P0finger t0Initial probability of collision at time, P1Is t0＜t≤tfA set of collision probabilities for the time interval. P0Expressed as the integral along the velocity on a spherical coordinate system:
Figure BDA0003075387070000091
pg(r0；μr(t0),Pr(t0) ) mean value of μr(t0) Covariance of Pr(t0) A multi-dimensional probability density distribution function. When two meeting targets are far away from each other P0＝0。
Figure BDA0003075387070000092
Wherein
Figure BDA0003075387070000093
σ2＝r’T(Pv-PvrPr -1Pvr T)r’
Figure BDA0003075387070000094
Where the spherical integration is performed by Lebedev integration, the integration over time may be performed using standard integration equations, such as Newton-Cotes, Gauss-Konrod and Clenshaw-Curtis.
Therefore, the collision probability calculation process of the Gaussian mixture covariance evolution is completed.
Example 2
The embodiment 2 of the invention provides a collision probability calculation system of Gaussian mixture covariance evolution, which specifically processes the method of the embodiment 1, and the system comprises the following steps: the system comprises a Gaussian mixture method orbit covariance fitting module, a covariance evolution module and a collision probability calculation module; wherein the content of the first and second substances,
the Gaussian mixture method orbit covariance fitting module is used for performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit;
the covariance evolution module is used for carrying out covariance evolution on each Gaussian unit to obtain corresponding orbit evolution error distribution expressed by Gaussian mixture;
and the collision probability calculation module is used for calculating collision probability among the Gaussian units according to the orbit evolution error distribution expressed by the Gaussian mixture, and performing weighted summation to obtain the total collision probability.
Example 3
A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the method of embodiment 1 when executing the computer program.
Example 4
A computer-readable storage medium storing a computer program which, when executed by a processor, causes the processor to perform the method of embodiment 1.
Simulation example
A case is introduced into simulation analysis to verify that the Gaussian mixture method is used for forecasting the covariance and calculating the collision probability, the initial covariance of the case is large, and the track state and the covariance are respectively forecasted by a linearization method, an unscented Kalman filtering method, a Gaussian mixture method + UT method and a Monte Carlo simulation method in the case. The forecast interval was 2 days. The Monte Carlo method simulates one million sample points containing random errors according to the initial orbit covariance, and each sample point is forecasted forwards for 2 days to obtain an orbit sample distribution graph at the forecasting time. The Monte Carlo method is the most accurate covariance forecast method and is used as a basis for verifying the precision of other methods. The covariance calculated by the forecasting simulation method is applied to the calculation of the collision probability of two targets, and the collision probability under different covariance forecasting methods is obtained.
FIG. 2 is a 3 σ error contour for radial along-the-track planar projection of Monte Carlo samples and linear, unscented Kalman Filter, and Gaussian mixture methods. After 2 days of initial orbit covariance evolution, the distribution of monte carlo sample points no longer obeys the gaussian distribution, but instead assumes a crescent-shaped distribution in the radial direction. The linearization method and the unscented kalman filter method still assume that the orbit state distribution is gaussian. The unscented kalman filter method can describe the overall distribution area of the orbit state, but cannot represent a specific distribution shape. The linearization method cannot even represent the distribution area of the track state. Therefore, the covariance of the linearization method and unscented kalman filter evolution cannot capture the orbit covariance distribution characteristic. And the Gaussian mixture distribution has better fitting effect with the distribution of Monte Carlo sample points. Table 1 shows the results of collision probability calculations for different methods. The collision probability calculated by the linearization method and the unscented kalman filter method is smaller than that calculated by the monte carlo method, and the calculation precision of the method cannot be used as a standard for measuring the space target collision early warning. The collision probability result of the Gaussian mixture calculation shows that the calculation accuracy of the collision probability is improved by one order of magnitude by the Gaussian mixture covariance calculation method.
TABLE 1 Collision probability calculation results for different methods
Figure BDA0003075387070000111
Finally, it should be noted that the above embodiments are only used for illustrating the technical solutions of the present invention and are not limited. Although the present invention has been described in detail with reference to the embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (8)

1. A method of collision probability computation of gaussian mixture covariance evolution, the method comprising:
performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit;
carrying out covariance evolution on each Gaussian unit to obtain corresponding orbital evolution error distribution represented by Gaussian mixture;
and calculating collision probability among Gaussian units aiming at the orbit evolution error distribution represented by the Gaussian mixture, and performing weighted summation to obtain the total collision probability.
2. The method for calculating collision probability of Gaussian mixture covariance evolution according to claim 1, wherein the two spatial objects to be predicted are subjected to initial orbit covariance fitting based on a Gaussian mixture algorithm respectively to obtain a Gaussian mixture probability density function of each spatial object with respect to a Gaussian unit; the method specifically comprises the following steps:
respectively calculating Gaussian probability density of each space object relative to the initial orbit based on the divided Gaussian units; fitting the covariance based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function P (x) of the space object, wherein the Gaussian mixture probability density function P (x) is as follows:
Figure FDA0003075387060000011
where N is the total number of Gaussian units, pgi(x；mi,Pi) Is the probability representation of the ith Gaussian unit, x represents the state quantity of the initial orbit and comprises a three-dimensional velocity vector and a three-dimensional position vector, miIs the mean value of the probability density of the ith Gaussian cell, PiIs the covariance of the probability density of the ith Gaussian cell, alphaiIs the weighted value of the probability density function of the ith Gaussian cell.
3. The method for computing collision probability of gaussian mixture covariance evolution according to claim 2, wherein the covariance evolution of each gaussian unit is performed to obtain an orbital evolution error distribution represented by a corresponding gaussian mixture; the method specifically comprises the following steps:
and carrying out covariance evolution on each Gaussian unit of each space object by adopting an Unscented transformation method, taking the estimated value and covariance of the state quantity x of the initial orbit at the 0 th moment as initial values, taking a discretization nonlinear equation as a model, and carrying out Unscented transformation until the estimated value and covariance at the T th moment are obtained, so that propagation of the mean value and covariance of the Gaussian units is realized to obtain the orbit evolution error distribution expressed by Gaussian mixture.
4. The method for calculating collision probability of Gaussian mixture covariance evolution according to claim 3, wherein the estimated value and covariance of the state quantity x of the initial orbit at the 0 th time are used as initial values, a discretized nonlinear equation is used as a model, and the Unscented transformation is performed until the estimated value and covariance at the T th time are obtained; the method specifically comprises the following steps:
according to the estimated value of the state quantity x of the ith Gaussian unit initial orbit at the k-1 moment
Figure FDA0003075387060000021
Sum covariance
Figure FDA0003075387060000022
1≤k≤T：
Figure FDA0003075387060000023
Figure FDA0003075387060000024
Build Sigma dots, get the following formula:
Figure FDA0003075387060000025
Figure FDA0003075387060000026
Figure FDA0003075387060000027
wherein σk-1Is and
Figure FDA0003075387060000028
a column vector of the same type is used,
Figure FDA0003075387060000029
the square root of the matrix is represented, lambda is an adjustable scale parameter and is used for improving approximation accuracy, and n is a state quantity dimension; sigma point χ is obtained by Cholesky decomposition or eigenvalue decomposition of covariance matrixk-1(j) Wherein Cholesky decomposition requires that the covariance matrix is positive, j represents the jth column of the covariance matrix;
directly carrying out state and observation prediction according to the calculated Sigma point according to the discretization nonlinear model to obtain a predicted Sigma sample with symmetric distributionDot chik(j) (j ═ 0, … 2n) is:
χk(j)＝f[χk-1(j),k]
wherein f [. cndot. ] represents a time-dependent nonlinear orbit prediction function, and k represents the kth moment;
according to Sigma sample point χk(j) Calculating mean of states
Figure FDA00030753870600000210
Sum variance
Figure FDA00030753870600000211
Figure FDA00030753870600000212
Figure FDA00030753870600000213
Wherein, χk(0) Represents the Sigma sample Point χk(j) In the 0 th column, alpha and beta are constant values, and the value range of alpha is more than or equal to 1 and less than or equal to 10-4In the case of the gaussian distribution, β is 2.
5. The method for calculating the collision probability of the Gaussian mixture covariance evolution as claimed in claim 4, wherein the collision probability is calculated for the orbital evolution error distribution represented by the Gaussian mixture through the Gaussian units, and the total collision probability is obtained by weighted summation; the method specifically comprises the following steps:
the time interval t is obtained according to the following formula0→tfProbability of collision of fcComprises the following steps:
fc＝P0+P1
wherein, P0Is t0Initial probability of collision at time, P1Is t0＜t≤tfA set of collision probabilities for the time intervals;
P0expressed as the integral of the velocity along the sphere coordinate system, satisfies the following equation:
Figure FDA0003075387060000031
wherein p isg(r0；μr(t0),Pr(t0) ) mean value of μr(t0) Covariance of Pr(t0) A multi-dimensional probability density distribution function, subscript r denotes the spherical radius, r0The radius of the spherical integral is represented, R represents the combined radius of the two intersecting objects, and theta and phi respectively represent the azimuth angle and the elevation angle of the spherical integral. When the distance between the space objects P and S exceeds a threshold value0＝0；
P1Satisfies the following formula:
Figure FDA0003075387060000032
where θ, φ represents the azimuth and elevation angles of the spherical integral, respectively, v (r ', t) represents the relative integral velocity, r' represents the integral radius, satisfying the following equation:
Figure FDA0003075387060000033
wherein the content of the first and second substances,
Figure FDA0003075387060000034
is an expression of the original variable, and satisfies the following formula:
Figure FDA0003075387060000035
Figure FDA0003075387060000036
σ2＝r’T(Pv-PvrPr -1Pvr T)r’
wherein, PvCovariance P representing the ellipsoid of the joint errorrelThe 4 th to 6 th column parts, P, corresponding to the 4 th to 6 th rowsvrCovariance P representing the ellipsoid of the joint errorrelColumn 1 to column 3 portions corresponding to row 4 to row 6; η represents an integral variable;
a Lebedev integration method is adopted for spherical integration, and the integration changing along with time adopts a standard integration formula;
the total collision probability P is obtained according to the following formulacComprises the following steps:
Figure FDA0003075387060000041
wherein f isc(. is the probability of collision between two Gaussian mixture units, NP,NSRespectively representing the number of Gaussian units of the space object P and the space object S, i and j respectively representing the ith Gaussian unit of the space object P and the jth Gaussian unit of the space object S,
Figure FDA0003075387060000042
respectively representing the weighting coefficient of the ith Gaussian unit of the spatial object P, the weighting coefficient of the jth Gaussian unit of the spatial object S, and TCA (i, j) representing the intersection time of the ith Gaussian unit and the jth Gaussian unit.
6. A collision probability calculation system of gaussian mixture covariance evolution, characterized in that the system comprises: the system comprises a Gaussian mixture method orbit covariance fitting module, a covariance evolution module and a collision probability calculation module; wherein the content of the first and second substances,
the Gaussian mixture method orbit covariance fitting module is used for performing initial orbit covariance fitting on two space objects to be predicted respectively based on a Gaussian mixture algorithm to obtain a Gaussian mixture probability density function of each space object relative to a Gaussian unit;
the covariance evolution module is used for carrying out covariance evolution on each Gaussian unit to obtain corresponding orbit evolution error distribution expressed by Gaussian mixture;
and the collision probability calculation module is used for calculating collision probability among the Gaussian units according to the orbit evolution error distribution expressed by the Gaussian mixture, and performing weighted summation to obtain the total collision probability.
7. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, characterized in that the processor implements the method of any of claims 1 to 5 when executing the computer program.
8. A computer-readable storage medium, characterized in that the computer-readable storage medium stores a computer program which, when executed by a processor, causes the processor to carry out the method of any one of claims 1 to 5.
CN202110550985.5A 2021-05-20 2021-05-20 Method and system for calculating collision probability of Gaussian mixture covariance evolution Active CN113343436B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110550985.5A CN113343436B (en) 2021-05-20 2021-05-20 Method and system for calculating collision probability of Gaussian mixture covariance evolution

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110550985.5A CN113343436B (en) 2021-05-20 2021-05-20 Method and system for calculating collision probability of Gaussian mixture covariance evolution

Publications (2)

Publication Number Publication Date
CN113343436A true CN113343436A (en) 2021-09-03
CN113343436B CN113343436B (en) 2022-02-18

Family

ID=77470153

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110550985.5A Active CN113343436B (en) 2021-05-20 2021-05-20 Method and system for calculating collision probability of Gaussian mixture covariance evolution

Country Status (1)

Country Link
CN (1) CN113343436B (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114329887A (en) * 2021-11-22 2022-04-12 北京理工大学 Data-driven dynamic Gaussian mixture spacecraft orbit error evolution method
CN114973783A (en) * 2022-08-02 2022-08-30 中国人民解放军63921部队 Spatial target collision early warning criterion optimization method

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050255253A1 (en) * 2004-05-13 2005-11-17 White John M Apparatus and methods for curing ink on a substrate using an electron beam
CN106249756A (en) * 2016-09-20 2016-12-21 北京理工大学 A kind of planetary landing obstacle avoidance control method based on collision probability
US20170124758A1 (en) * 2015-10-29 2017-05-04 Intelligent Fusion Technology, Inc Method and system for predicting collision probability of space objects via graphics processing unit
CN107544067A (en) * 2017-07-06 2018-01-05 西北工业大学 One kind is based on the approximate Hypersonic Reentry Vehicles tracking of Gaussian Mixture
CN108279011A (en) * 2018-01-30 2018-07-13 北京理工大学 Planetary detection landing path comprehensive optimization method
CN108507593A (en) * 2018-04-09 2018-09-07 郑州轻工业学院 A kind of dimensionality reduction RTS ellipsoid collection person's smoothing methods of INS errors model
CN110379203A (en) * 2019-06-21 2019-10-25 江苏大学 A kind of driving steering anti-collision warning method
CN111428369A (en) * 2020-03-26 2020-07-17 中国人民解放军32035部队 Method for calculating confidence of space target collision early warning result

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050255253A1 (en) * 2004-05-13 2005-11-17 White John M Apparatus and methods for curing ink on a substrate using an electron beam
US20170124758A1 (en) * 2015-10-29 2017-05-04 Intelligent Fusion Technology, Inc Method and system for predicting collision probability of space objects via graphics processing unit
CN106249756A (en) * 2016-09-20 2016-12-21 北京理工大学 A kind of planetary landing obstacle avoidance control method based on collision probability
CN107544067A (en) * 2017-07-06 2018-01-05 西北工业大学 One kind is based on the approximate Hypersonic Reentry Vehicles tracking of Gaussian Mixture
CN108279011A (en) * 2018-01-30 2018-07-13 北京理工大学 Planetary detection landing path comprehensive optimization method
CN108507593A (en) * 2018-04-09 2018-09-07 郑州轻工业学院 A kind of dimensionality reduction RTS ellipsoid collection person's smoothing methods of INS errors model
CN110379203A (en) * 2019-06-21 2019-10-25 江苏大学 A kind of driving steering anti-collision warning method
CN111428369A (en) * 2020-03-26 2020-07-17 中国人民解放军32035部队 Method for calculating confidence of space target collision early warning result

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
KYLEJ.DEMARS 等: "Collision Probability with Gaussian Mixture Orbit Uncertainty", 《JOURNAL OF GUIDANCE，CONTROL，AND DYNAMICS》 *
VIVEK VITTALDEV 等: "Collision probability for space objects using Gaussian mixture models", 《ADVANCES IN THE ASTRONAUTICAL SCIENCES》 *
张莹 等: "基于协方差分析法的轨控误差对碰撞概率影响研究", 《第十一届中国卫星导航年会论文集——S04 卫星轨道与系统误差处理》 *
白显宗: "空间目标轨道预报误差与碰撞概率问题研究", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 *
闫瑞东 等: "基于数值轨道模型的轨道协方差演化分析", 《红外与激光工程》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114329887A (en) * 2021-11-22 2022-04-12 北京理工大学 Data-driven dynamic Gaussian mixture spacecraft orbit error evolution method
CN114973783A (en) * 2022-08-02 2022-08-30 中国人民解放军63921部队 Spatial target collision early warning criterion optimization method

Also Published As

Publication number Publication date
CN113343436B (en) 2022-02-18

Similar Documents

Publication Publication Date Title
CN111985093B (en) Adaptive unscented Kalman filtering state estimation method with noise estimator
CN106054170B (en) A kind of maneuvering target tracking method under constraints
Gultekin et al. Nonlinear Kalman filtering with divergence minimization
CN111722214B (en) Method for realizing radar multi-target tracking PHD
CN113343436B (en) Method and system for calculating collision probability of Gaussian mixture covariance evolution
Li et al. Geometry-driven deterministic sampling for nonlinear bingham filtering
Sun et al. Modeling and tracking of maneuvering extended object with random hypersurface
CN112800599B (en) Non-grid DOA estimation method based on ADMM under array element mismatch condition
CN113221311A (en) Uncertainty quantification method for wind speed of atmospheric boundary layer
CN116047495B (en) State transformation fusion filtering tracking method for three-coordinate radar
CN111708013B (en) Target tracking filtering method for distance coordinate system
CN117036400A (en) Multi-target group tracking method based on fuzzy clustering data association of Gaussian mixture model
Xiao et al. A multiple model particle filter for maneuvering target tracking based on composite sampling
CN116224320B (en) Radar target tracking method for processing Doppler measurement under polar coordinate system
Jia et al. Adaptive cubature Kalman filter with directional uncertainties [Correspondence]
CN115544425A (en) Robust multi-target tracking method based on target signal-to-noise ratio characteristic estimation
CN111722213B (en) Pure distance extraction method for maneuvering target motion parameters
Godsill Particle filters for continuous-time jump models in tracking applications
CN113065287B (en) Small celestial body gravitational field rapid prediction method based on implicit characteristics
Barlow et al. An alternative algorithm for the refinement of ULV decompositions
Chen et al. Improved Unscented Kalman Filtering Algorithm Applied to On-vehicle Tracking System
Wang et al. A moving grid cell based MCL algorithm for mobile robot localization
CN113688498B (en) Method for determining trusted interval response of aircraft structure
Zhuang et al. Accurate statistical modeling method for dynamic RCs
CN112115419B (en) System state estimation method and system state estimation device

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant