Space target optical segmental arc rapid association method based on point cloud intersection

Abstract

The invention discloses a point cloud intersection-based space target optical segmental arc rapid association method. The method comprises the following steps: 1, acquiring optical segmental arc data and constructing a constraint admissible domain; 2, generating a track point cloud; 3, forecasting parallel orbit point cloud; 4, fast point cloud intersection is carried out; and 5, determining an initial orbit based on point cloud intersection. Through the above process, the space target rapid segmental arc correlation method based on point cloud intersection is provided, an admissible domain is constructed by two optical observation segmental arcs and orbit point clouds are generated, and then the point clouds of the two segmental arcs are forecasted to the same moment through a parallel orbit point cloud forecasting technology. Quickly locking an overlapping region of the two point clouds by adopting a minimum bounding box method, and judging the relevance between the point clouds through a mahalanobis distance; and finally, based on the Bayesian estimation principle, carrying out initial orbit determination on the successfully associated arc segments.

Description

Technical Field

[0001] The present invention provides a method for quickly associating optical arc segments of space targets based on point cloud intersection, which involves using an admissible region method and a minimum bounding box method to judge the correlation between two optical observation arc segments and belongs to the field of space situation awareness. Background Art

[0002] With the rapid development of the aerospace industry, particularly the advancement of mega-constellations, the number of space objects in orbit has exploded in recent years. By the end of 2024, the US space detection system had cataloged 29,000 space objects. Limited by current monitoring technology, the existing cataloging system primarily targets space objects larger than 10 centimeters. However, even space debris with a diameter of only 1 centimeter can have enough impact energy to cause fatal damage to operating satellites. The European Space Debris Office estimates that as of August 1, 2024, there will be approximately 1.2 million space objects between 1 and 10 centimeters in size in space. Therefore, a large number of uncataloged small space debris poses a significant threat to satellites in orbit.

[0003] Space targets are typically tracked using ground-based radar, ground-based telescopes, and space-based telescopes. Space-based telescopes are unaffected by day and night, as well as atmospheric interference, and can significantly enhance the spatial and temporal coverage of space surveillance. Consequently, the use of space-based telescopes to detect small debris has attracted considerable attention in recent years. When space debris passes through the telescope's field of view, the system acquires a series of angular measurements of right ascension, declination, and other aspects through continuous imaging and astronomical positioning algorithms. These continuous observations constitute an optical arc segment. Space-based optical arc segments are typically short-arc observations. Due to the lack of distance information, the orbit determination accuracy of a single short arc segment is limited. Therefore, it is necessary to correlate multiple arc segments belonging to the same space target.

[0004] In early studies of optical arc correlation, correlation methods based on admissible domains faced significant efficiency bottlenecks and accuracy limitations. Traditional methods typically construct a Delaunay triangulation within the admissible domain and use the Covariance-Based Tracklet Correlation (CBTA) method to perform correlation tests at each grid point. This approach has two major drawbacks: First, to ensure that valid information is not lost, the grid must be densely divided. With optical monitoring systems generating tens of thousands of observation arcs daily, the computational complexity of this exhaustive search becomes prohibitive. Second, while some studies have attempted to implement admissible domain intersection strategies within the Delaunay state space, these strategies still rely on manual judgment of intersection and lack a systematic, efficient intersection algorithm. Other studies have mapped the admissible domain to the Poincaré orbital element space and evaluated the intersection of the admissible domains using probability density distributions. However, this approach requires computing integrals in the probability density space, which is computationally complex. Crucially, most existing studies only consider the average effects of the Earth's gravitational field and the J2 perturbation. In reality, orbital perturbations in the low-Earth orbit (LEO) environment (such as atmospheric drag and the gravitational pull of the Sun and Moon) have significant influence. This makes the application of intersection methods based on simplified dynamic models in LEO difficult to assess. Therefore, developing an arc-segment correlation algorithm that can adapt to complex perturbation environments while also offering efficient computational performance has become a pressing technical challenge.

[0005] In summary, the present invention proposes a point cloud intersection arc segment association method based on parallel computing and intelligent screening. A probabilistic orbit point cloud is constructed through a Gaussian mixture model, and the traditional deterministic association is transformed into a probabilistic matching. A parallel integration algorithm based on the improved vernal equinox orbit element is developed to achieve rapid propagation of large-scale orbit point clouds. At the same time, a minimum bounding box screening strategy is designed to significantly improve the efficiency of point cloud intersection. An accurate perturbation model and an adaptive covariance propagation mechanism are introduced to ensure the accuracy of association in complex low-orbit environments. This method breaks through the efficiency bottleneck of traditional arc segment association technology and provides a reliable technical means for real-time monitoring of massive space targets. Summary of the Invention

[0006] (1) Purpose of the invention

[0007] This invention aims to provide a method for rapidly associating optical arc segments of space targets based on point cloud intersection. By constructing a probabilistic orbit point cloud model and combining parallel computing with intelligent screening techniques, this method achieves efficient association of optical observation data of space targets. This method primarily addresses the problems of low computational efficiency and insufficient association accuracy encountered in existing technologies when processing large-scale orbit data. By leveraging innovative techniques such as Gaussian mixture models to characterize orbital uncertainty, improved parallel integration of equinox orbital elements, and multi-level minimum bounding box screening, this method significantly improves the speed and accuracy of arc segment association, providing reliable technical support for space target cataloging and collision warning.

[0008] (2) Technical solution

[0009] The present invention relates to a method for rapidly associating optical arcs of space targets based on point cloud intersection. This method requires two optical observation arcs observed by a low-orbit observation satellite or a ground-based optical observation station as input data to determine the correlation between the two arcs. First, the corresponding right ascension and declination angles and angle change rates of the two observation arcs are extracted, and a distance-distance change rate tolerance domain is constructed to generate an orbit point cloud. Then, the orbit point clouds of the two arcs are predicted to the same time using parallel orbit point cloud prediction technology. Next, the minimum bounding box method is used to quickly locate the overlapping area of ​​the two point clouds, and the correlation between the point clouds is determined using the Mahalanobis distance. Finally, the initial orbit of the successfully associated arc is determined based on the Bayesian estimation principle.

[0010] The method for quickly associating optical arc segments of space targets based on point cloud intersection described in the present invention is as follows: Figure 1 , and its implementation steps are as follows:

[0011] Step 1: Obtain optical arc segment data and construct the constraint allowable domain

[0012] Given an optical arc segment, the admissible region method converts continuous angle measurements into angle and angular rate measurements at intermediate moments in the arc segment:

[0013]

[0014] Generally called is the arc segment characterization vector, α, β are the measured right ascension and declination, is the rate of change of right ascension and declination angles. Figure 2 The observation model of , the sight vector is:

[0015]

[0016] The distance between the space target and the low-orbit monitoring platform is ρ, and the distance change rate is It is known that the position and velocity of the optical measuring device at this moment are r o With v0, given a set of distances and distance change rates Then the position of the space target r t With speed v t for

[0017]

[0018] in,

[0019]

[0020] So the orbital state can be expressed as For two arc segments observed at time t1 and time t2, respectively, given a set of At this point, the states of the first and last orbits are known. Space targets that are in low orbit and can be observed must meet the constraints of observation distance, semi-major axis, eccentricity, perigee height, etc., so it is possible to establish and The two-dimensional constraint allowed domain of , each value point in the allowed domain represents an orbit.

[0021] Step 2: Track point cloud generation

[0022] From the perspective of probability, the probability of any point in the admissible domain being a true value is the same, so it is believed that The joint probability density of The interior is uniformly distributed, which can be described mathematically as

[0023]

[0024] Where S is the allowed domain If the area is approximated by a discretized grid, then

[0025]

[0026] Where N is the number of grid points, Δρ and The sampling intervals of distance and distance change rate are respectively. Using Gaussian mixture model to approximate the uniform distribution characteristics inside the admissible domain, any grid point inside the admissible domain is regarded as the center of a two-dimensional Gaussian sample, then

[0027]

[0028] in, Represents normal distribution, and the calculation expression is

[0029]

[0030] Where m and P represent the expectation and covariance of the orbital state, respectively. Combined with the angle / angle change rate information, any point in the admissible domain represents an orbit, which is described by the state vector

[0031]

[0032] A point in the orbit point cloud corresponding to the discretized allowed domain is χ i The covariance of

[0033]

[0034] Step 3: Parallel track point cloud prediction

[0035] Compared with the Kepler orbital elements [a,e,i,ω,Ω,θ], the improved equinox orbital elements ξ=[p,f,g,h,j,k,L] T There is no singularity at zero eccentricity, defined as

[0036]

[0037] Where a, e, i, ω, Ω, and θ are the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and true anomaly, respectively, as defined by Kepler orbital elements. For prograde orbits (i ≤ 90°), κ = 1; for retrograde orbits (i > 90°), κ = -1, and L is called true longitude. The improved Lagrangian planetary equations for the vernal equinox orbital elements are:

[0038]

[0039] Among them, U represents the perturbation potential function, s 2 and w are intermediate variables, and the expression is

[0040] s 2 =1+h 2 +k 2 (13)

[0041] w=1+fcosL+gsinL (14)

[0042] The improved differential equation for the evolution of the vernal equinox orbital elements is obtained

[0043]

[0044] Given the initial value, the ordinary differential equation can be used to predict the state by numerical integration. If the evolution of multiple orbits is described at the same time, it can be written as

[0045]

[0046] in,

[0047] [ξ]=[ξ1,ξ2,...,ξ n ](17)

[0048]

[0049] Among them, ξ i and the χ obtained by the admissible region method i One-to-one correspondence. Direct calculation using Hadamard product The Hadamard product is also called the dot multiplication operation. The Hadamard product of two matrices or vectors of the same dimension is

[0050]

[0051] Calculate integrals using Hadamard products Parallel track point cloud prediction can be achieved.

[0052] In addition, in order to reflect the statistical distribution of the uncertainty of the initial orbit determination, each sample point in the orbit point cloud has an uncertainty, so the uncertainty of each orbit also needs to be predicted. Taking a single orbit as an example, using linear mapping approximation, the disturbance at the initial time χ0 causes the final time ξ f The change is

[0053] δξ f =Φ χ2ξ (t f ,t0)δχ0(20)

[0054] Among them, Φ χ2ξ (t f ,t0) is the mapping matrix. The uncertainty of χ0 at the initial moment is P χ , then at the final moment ξ f The uncertainty is

[0055] P ξ =Φ χ2ξ (t f ,t0)P χ [Φ χ2ξ (t f ,t0)] T (twenty one)

[0056] Therefore, the core of orbit uncertainty lies in calculating Φ χ2ξ (t f ,t0). Parallel prediction of track point cloud state values ​​has been achieved, and the Φ of different tracks χ2ξ (t f ,t0) is difficult to implement. The present invention directly uses the numerical difference method to calculate Φ χ2ξ (t f ,t0). The track point cloud is predicted seven times in parallel. The first time is to obtain the state prediction value, and the remaining six times are to obtain Φ χ2ξ (t f ,t0), and then the predicted orbit uncertainty is obtained.

[0057] Step 4: Quick point cloud intersection

[0058] If two arcs belong to the same space target, then the two admissible domains obtained contain the same orbit. Without considering observation errors and model errors, the two admissible domains are predicted to the same time, and the two orbital state sets should intersect at the true value of the orbit. And because the orbit has only six degrees of freedom and the two admissible domains provide eight constraints, it is mathematically expressed as the solution of an overdetermined equation. In general, there is a unique solution or no solution, and there is only one intersection point corresponding to the two orbital state sets.

[0059] In order to realize the fast point cloud intersection operation, the present invention adopts the Figure 4 The minimum bounding box method is used to quickly locate the overlapping area of ​​the two point clouds. The two point clouds formed by the two arc segments are called C1 and C2. The physical quantity X of each data point in C1 is extracted, such as the semi-path p, the improved vernal equinox orbit elements f, g, h, k, and the maximum and minimum values ​​are X respectively. 1max 、X 1min Similarly, extract the X of each data point in C2 and get its maximum and minimum values ​​X 2max 、X 2min , the maximum and minimum values ​​of the minimum bounding box are min{X 1max ,X 2max}、max{X 1min ,X 2min}, obviously, if min{X 1max ,X 2max}<max{X 1min ,X 2min}, then there is no minimum bounding box and the two point clouds do not intersect. If there is a minimum bounding box, discard the data points outside the minimum bounding box in C1 and C2, and you will get X 1min =X 2min 、X 1max =X 2max .

[0060] After multiple rounds of iterations to find the minimum bounding box for p, f, g, h, k, the number of data points will be significantly reduced, and even converge to a very small area containing the intersection of the two point clouds. Taking into account the observation error, the dynamic model error, and the sampling interval of the distance / distance change rate, the minimum bounding box can be appropriately enlarged. For example, the minimum and maximum values ​​of the minimum bounding box can be taken as min{X 1max ,X 2max}+ΔX、max{X 1min ,X 2min In addition, the present invention does not perform the minimum bounding box operation on the true longitude L, because the true longitude L is a fast variable and is not as stable as other orbital elements.

[0061] The result of the initial orbit determination in the admissible domain is described by the orbit point cloud and the uncertainty of each point, that is, the Gaussian mixture model. The true orbit value falls within at least one Gaussian sample. If two arcs belong to the same spatial target, there is a pair of Gaussian samples that "intersect". Whether a pair of Gaussian samples "intersect" can be statistically judged using the Mahalanobis distance, which is defined as

[0062]

[0063] Among them, X1 and P1 are the mean and covariance of the first Gaussian sample, and X2 and P2 are the mean and covariance of the second Gaussian sample. If two Gaussian samples intersect, χ 2 Should satisfy the chi-square distribution with 6 degrees of freedom, when χ 2 <16.82 with a confidence level of 99%. In other words, if there is a pair of Gaussian samples with a Mahalanobis distance less than It is considered that the two arcs may belong to the same spatial target.

[0064] After obtaining the overlapping area of ​​two track point clouds using the minimum bounding box operation, each point cloud only contains a small number of sample points. Figure 5 The remaining sample points are combined in pairs and the Mahalanobis distances of all combinations are calculated. If there is a pair of combinations with a Mahalanobis distance less than 4.1, it is considered that the two arcs may belong to the same spatial target.

[0065] Step 5: Determine the initial trajectory based on point cloud intersection

[0066] If the two arc segments obtained after the point cloud intersection are related, the two point clouds after the minimum bounding box operation are written as Gaussian mixture models

[0067]

[0068] in, and is the mean and covariance of the i-th Gaussian sample in the first point cloud, and is the mean and covariance of the i-th Gaussian sample in the second point cloud. Taking the first point cloud as prior information and the second point cloud as observation, according to the Bayesian estimation principle, we get

[0069]

[0070] in

[0071]

[0072] Further

[0073]

[0074] in

[0075]

[0076]

[0077] Finally, Equation (25) can be written as

[0078]

[0079] Among them, w ij is the weight coefficient:

[0080]

[0081] The initial orbit determination results of the two arc segments are

[0082]

[0083] In summary, the implementation steps of the present invention are as follows: Figure 1 As shown in the figure. The optical arc segment association and preliminary orbit determination method for low-orbit space targets proposed through the above process constructs a new orbit point cloud processing framework based on admissible domain theory and Gaussian mixture models, which is fundamentally different from traditional methods. Its core technical concept is to transform the uncertainty of optical observations into orbit point clouds by constructing a constrained admissible domain, utilize parallel computing technology to achieve efficient propagation of large-scale orbit point clouds, and innovatively combine the minimum bounding box method with the Mahalanobis distance to achieve fast point cloud matching. Finally, the preliminary orbit is obtained by fusing the associated arc segment information through Bayesian estimation.

[0084] (3) Advantages

[0085] The advantages of the method for quickly associating optical arc segments of space targets based on point cloud intersection provided by the present invention are:

[0086] ① The optical arc segment correlation method proposed in the present invention is applicable to free-moving space targets in various orbital environments and has a wide range of applications;

[0087] ② The present invention adopts Gaussian mixture model to describe the orbital distribution in the admissible domain, which can maintain the continuous probability characteristics.

[0088] It is also convenient for numerical calculations and provides a mathematical basis for subsequent statistical analysis;

[0089] ③ Based on the improved parallel propagation algorithm of the vernal equinox orbit elements, the present invention uses the Hadamard product to achieve efficient prediction of tens of thousands of orbit point clouds;

[0090] ④ This invention innovatively combines the minimum bounding box iterative shrinkage method and the Mahalanobis distance test, which can quickly lock potential associated pairs in massive track point clouds, greatly reducing the computational complexity and improving the association efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 It is a flow chart of the implementation steps of the present invention.

[0092] Figure 2 This is a schematic diagram of the observation model in this method.

[0093] Figure 3 It is a schematic diagram of the constraint allowed domain.

[0094] Figure 4 It is a diagram of the minimum bounding box.

[0095] Figure 5 This is a schematic diagram of the point cloud distribution of two tracks after the minimum bounding box operation. DETAILED DESCRIPTION

[0096] The specific implementation process of the present invention will be further described in detail below in conjunction with the technical solution.

[0097] The present invention relates to a method for rapidly associating optical arcs of space targets based on point cloud intersection. This method requires two optical observation arcs observed by a low-orbit observation satellite or a ground-based optical observation station as input data to determine the correlation between the two arcs. First, the corresponding right ascension and declination angles and angle change rates of the two observation arcs are extracted, and a distance-distance change rate tolerance domain is constructed to generate an orbit point cloud. Then, the orbit point clouds of the two arcs are predicted to the same time using parallel orbit point cloud prediction technology. Next, the minimum bounding box method is used to quickly locate the overlapping area of ​​the two point clouds, and the correlation between the point clouds is determined using the Mahalanobis distance. Finally, the initial orbit of the successfully associated arc is determined based on the Bayesian estimation principle.

[0098] The method for quickly associating optical arc segments of space targets based on point cloud intersection described in the present invention is as follows: Figure 1 , and its implementation steps are as follows:

[0099] Step 1: Obtain optical arc segment data and construct the constraint allowable domain

[0100] Given an optical arc segment, the admissible region method converts continuous angle measurements into angle and angular rate measurements at intermediate moments in the arc segment:

[0101]

[0102] Generally called is the arc segment characterization vector, α, β are the measured right ascension and declination, is the rate of change of the angle of right ascension and declination. The calculation method is usually quadratic curve fitting, assuming that the change of angle with time is a quadratic curve relationship, that is,

[0103]

[0104]

[0105] make

[0106]

[0107] Assuming that the measurement noise level at each moment in each time period is similar, using the least squares method, we can get

[0108]

[0109] And the covariance matrix of the extracted angle / angular rate can be obtained, which is recorded as

[0110] Refer to the attached Figure 2 The observation model of , the sight vector is:

[0111]

[0112] The distance between the space target and the low-orbit monitoring platform is ρ, and the distance change rate is It is known that the position and velocity of the optical measuring device at this moment are r o With v0, given a set of distances and distance change rates Then the position of the space target r t With speed v t for

[0113]

[0114] in

[0115]

[0116] So the orbital state can be expressed as For two arc segments observed at time t1 and time t2, respectively, given a set of At this point, the states of the first and last orbits are known. Space targets that are in low orbit and can be observed must meet the constraints of observation distance, semi-major axis, eccentricity, perigee height, etc., so it is possible to establish and The two-dimensional constraint allowed domain is , and each value point in the allowed domain represents an orbit. Take a LEO target with an orbital altitude of 1000km as an example, refer to the attached Figure 3 To illustrate the constraints of the allowed domain, the sampling intervals of distance and distance change rate are 20 km and 20 m / s respectively. The constraints are: energy constraint ∈ < 0, semi-major axis constraint 6678 km ≤ a ≤ 8400 km, eccentricity constraint e ≤ 0.2, perigee constraint rp >6678km.

[0117] Step 2: Track point cloud generation

[0118] From the perspective of probability, the probability of any point in the admissible domain being a true value is the same, so it is believed that The joint probability density of The interior is uniformly distributed, which can be described mathematically as

[0119]

[0120] Where S is the allowed domain If the area is approximated by a discretized grid, then

[0121]

[0122] Where N is the number of grid points, Δρ and The sampling intervals of distance and distance change rate are respectively. Using Gaussian mixture model to approximate the uniform distribution characteristics inside the admissible domain, any grid point inside the admissible domain is regarded as the center of a two-dimensional Gaussian sample, then

[0123]

[0124] in, Represents normal distribution, and the calculation expression is

[0125]

[0126] Where m and P represent the expectation and covariance of the orbital state, respectively. Combined with the angle / angle change rate information, any point in the admissible domain represents an orbit, which is described by the state vector

[0127]

[0128] A point in the orbit point cloud corresponding to the discretized allowed domain is χ i The covariance of

[0129]

[0130] Step 3: Parallel track point cloud prediction

[0131] Compared with the Kepler orbital elements [a,e,i,ω,Ω,θ], the improved equinox orbital elements ξ=[p,f,g,h,j,k,L] T There is no singularity at zero eccentricity, defined as

[0132]

[0133] Where a, e, i, ω, Ω, and θ are the semi-major axis, eccentricity, orbital inclination, right ascension of the ascending node, argument of perigee, and true anomaly, respectively, as defined by Kepler orbital elements. For prograde orbits (i ≤ 90°), κ = 1; for retrograde orbits (i > 90°), κ = -1, and L is called true longitude. The improved Lagrangian planetary equations for the vernal equinox orbital elements are:

[0134]

[0135] Among them, U represents the perturbation potential function, s 2 and w are intermediate variables, and the expression is

[0136] s 2 =1+h 2 +k 2 (50)

[0137] w=1+fcosL+gsinL (51)

[0138] If only the harmonic terms in the Earth's non-gravitational field are considered, the perturbation potential function is

[0139]

[0140] Where r is the distance from the center of the Earth, is the geocentric dimension, R e =6378136.6m is the radius of the earth, P n is the Legendre polynomial, J n To make the harmonic term sparse, we use the improved vernal equinox orbit element to get

[0141]

[0142] Further

[0143]

[0144] For LEO space targets, the four perturbations J2, J3, J4, and J6 have a more significant impact. In view of the requirement of fast arc segment association, this paper only considers these four perturbations. The corresponding Legendre polynomials are

[0145]

[0146] The improved differential equation for the evolution of the vernal equinox orbital elements is obtained

[0147]

[0148] Given the initial value, the ordinary differential equation can be used to predict the state by numerical integration. If the evolution of multiple orbits is described at the same time, it can be written as

[0149]

[0150] in

[0151] [ξ]=[ξ1,ξ2,...,ξ n ](62)

[0152]

[0153] Among them, ξ i and the χ obtained by the admissible region method i One-to-one correspondence. Direct calculation using Hadamard product The Hadamard product is also called the dot multiplication operation. The Hadamard product of two matrices or vectors of the same dimension is

[0154]

[0155] Calculate integrals using Hadamard products Parallel track point cloud prediction can be achieved.

[0156] In addition, in order to reflect the statistical distribution of the uncertainty of the initial orbit determination, each sample point in the orbit point cloud has an uncertainty, so the uncertainty of each orbit also needs to be predicted. Taking a single orbit as an example, using linear mapping approximation, the disturbance at the initial time χ0 causes the final time ξ f The change is

[0157] δξ f =Φ χ2ξ (t f ,t0)δχ0(65)

[0158] Among them, Φ χ2ξ (t f ,t0) is the mapping matrix. The uncertainty of χ0 at the initial moment is P χ , then at the final moment ξ f The uncertainty is

[0159] P ξ =Φ χ2ξ (t f ,t0)P χ [Φ χ2ξ (t f ,t0)] T (66)

[0160] Therefore, the core of orbit uncertainty lies in calculating Φ χ2ξ (t f ,t0). Parallel prediction of track point cloud state values ​​has been achieved, and the Φ of different tracks χ2ξ (t f,t0) is difficult to implement. The present invention directly uses the numerical difference method to calculate Φ χ2ξ (t f ,t0), generating six sets of orthogonal δχ0, corresponding to the right ascension / declination changes of 10 -6 rad, right ascension and declination rate of change 10 -8 rad / s, distance change 1km, distance change rate change 1m / s, after six orbit predictions, six sets of δξ are obtained f , Φ χ2ξ (t f ,t0) is calculated according to the following formula

[0161]

[0162] Therefore, the track point cloud is predicted seven times in parallel. The first time is to obtain the state prediction value, and the remaining six times are to obtain Φ χ2ξ (t f ,t0), and then the predicted orbit uncertainty is obtained.

[0163] Step 4: Quick point cloud intersection

[0164] If two arcs belong to the same space target, then the two admissible domains obtained contain the same orbit. Without considering observation errors and model errors, the two admissible domains are predicted to the same time, and the two orbital state sets should intersect at the true value of the orbit. And because the orbit has only six degrees of freedom and the two admissible domains provide eight constraints, it is mathematically expressed as the solution of an overdetermined equation. In general, there is a unique solution or no solution, and there is only one intersection point corresponding to the two orbital state sets.

[0165] In order to realize the fast point cloud intersection operation, the present invention adopts the Figure 4 The minimum bounding box method is used to quickly locate the overlapping area of ​​the two point clouds. The two point clouds formed by the two arc segments are called C1 and C2. The physical quantity X of each data point in C1 is extracted, such as the semi-path p, the improved vernal equinox orbit elements f, g, h, k, and the maximum and minimum values ​​are X respectively. 1max 、X 1min Similarly, extract the X of each data point in C2 and get its maximum and minimum values ​​X 2max 、X 2min , the maximum and minimum values ​​of the minimum bounding box are min{X 1max ,X 2max}、max{X 1min ,X 2min}, obviously, if min{X 1max ,X 2max}<max{X 1min ,X 2min}, then there is no minimum bounding box and the two point clouds do not intersect. If there is a minimum bounding box, discard the data points outside the minimum bounding box in C1 and C2, and you will get X 1min =X 2min 、X 1max =X 2max .

[0166] After multiple rounds of iterations to find the minimum bounding box for p, f, g, h, k, the number of data points will be significantly reduced, and even converge to a very small area containing the intersection of the two point clouds. Taking into account the observation error, the dynamic model error, and the sampling interval of the distance / distance change rate, the minimum bounding box can be appropriately enlarged. For example, the minimum and maximum values ​​of the minimum bounding box can be taken as min{X 1max ,X 2max}+ΔX、max{X 1min ,X 2min In addition, the present invention does not perform the minimum bounding box operation on the true longitude L, because the true longitude L is a fast variable and is not as stable as other orbital elements.

[0167] The result of the initial orbit determination in the admissible domain is described by the orbit point cloud and the uncertainty of each point, that is, the Gaussian mixture model. The true orbit value falls within at least one Gaussian sample. If two arcs belong to the same spatial target, there is a pair of Gaussian samples that "intersect". Whether a pair of Gaussian samples "intersect" can be statistically judged using the Mahalanobis distance, which is defined as

[0168]

[0169] Among them, X1 and P1 are the mean and covariance of the first Gaussian sample, and X2 and P2 are the mean and covariance of the second Gaussian sample. If two Gaussian samples intersect, χ 2 Should satisfy the chi-square distribution with 6 degrees of freedom, when χ 2 <16.82 with a confidence level of 99%. In other words, if there is a pair of Gaussian samples with a Mahalanobis distance less than It is considered that the two arcs may belong to the same spatial target.

[0170] After obtaining the overlapping area of ​​two track point clouds using the minimum bounding box operation, each point cloud only contains a small number of sample points. Figure 5 The remaining sample points are combined in pairs and the Mahalanobis distances of all combinations are calculated. If there is a pair of combinations with a Mahalanobis distance less than 4.1, it is considered that the two arcs may belong to the same spatial target.

[0171] Step 5: Determine the initial trajectory based on point cloud intersection

[0172] If the two arc segments obtained after the point cloud intersection are related, the two point clouds after the minimum bounding box operation are written as Gaussian mixture models

[0173]

[0174] in, and is the mean and covariance of the i-th Gaussian sample in the first point cloud, and is the mean and covariance of the i-th Gaussian sample in the second point cloud. Taking the first point cloud as prior information and the second point cloud as observation, according to the Bayesian estimation principle, we get

[0175]

[0176] in,

[0177]

[0178] Further

[0179]

[0180] in

[0181]

[0182] Finally, Equation (25) can be written as

[0183]

[0184] Among them, w ij is the weight coefficient:

[0185]

[0186] The initial orbit determination results of the two arc segments are

[0187]

Claims

1. A method for fast associating optical arc segments of space targets based on point cloud intersection, characterized by: The steps are as follows Step 1: Obtain optical arc segment data and construct the constraint allowable domain For a segment of optical observation arc data, extract the normalized vector of angle and angle change rate: Where α and β represent right ascension and declination; Represents the rate of change of the angle of right ascension and declination. The sight vector is: Let the coordinate vector of the observation satellite be r o , the velocity vector is v o , so the position vector and velocity vector of the target satellite are: where ρ is the relative distance, is the relative distance change rate, and: The orbital state can then be described as For two arc segments observed at time t1 and time t2, a set of At this point, the states of the first and last orbits are known. Space targets that are in low orbit and can be observed must meet the constraints of observation distance, semi-major axis, eccentricity, perigee height, etc., so it is possible to establish and The two-dimensional constraint allowed domain of , each value point in the allowed domain represents an orbit. Step 2: Track point cloud generation Discretize the admissible domain into grid points, and use the Gaussian mixture model to approximate the uniform distribution characteristics inside the admissible domain. Take any grid point inside the admissible domain as the center of a two-dimensional Gaussian sample, then we have Combined with the angle / angle change rate information, any point in the admissible domain represents a trajectory, which is described by the state vector A point in the orbit point cloud corresponding to the discretized allowed domain is χ i The covariance of Step 3: Parallel track prediction The improved orbital equinox elements {p,f,g,h,k,L} are used for dynamic modeling, J2, J3, J4, and J6 perturbations are considered, and batch orbital state parallel integration is achieved through Hadamard product. Combining one orbit prediction and six disturbance propagations, the propagation of state transfer matrix and covariance can be solved efficiently. Step 4: Fast point cloud matching For the variable X, X represents the semi-path p and the p, f, g, h, k in the improved equinox element, the minimum bounding box is iteratively screened, shrinking the bounding box boundary to min{X 1max ,X 2max }+ΔX、max{X 1min ,X 2min }-ΔX, and finally verify the relevance of the point cloud within the bounding box by the Mahalanobis distance criterion that meets the conditions: Step 5: Determine the initial trajectory based on point cloud intersection Based on the Bayesian principle, the initial trajectory is determined using two point clouds that have undergone minimum bounding box operations. Among them, w ij is the weight coefficient:

2. The method for rapid association of optical arc segments of space targets based on point cloud intersection according to claim 1, characterized in that: The specific implementation method of step 3 parallel orbit prediction is: Given the initial value, the ordinary differential equation can be used to predict the state by numerical integration. If the evolution of multiple orbits is described at the same time, it can be written as in [ξ]=[ξ1,ξ2,...,ξ n ] (14) Among them, ξ i and the χ obtained by the admissible region method i One-to-one correspondence. Use Hadamard product to calculate integral calculation Parallel track point cloud prediction can be achieved. In addition, in order to predict the uncertainty of each track, taking a single track as an example, the linear mapping approximation is adopted, and the disturbance at the initial time χ0 causes the disturbance at the final time ξ f The change is right f =Φ χ2ξ (t f ,t0)δχ0 (16) Among them, Φ χ2ξ (t f ,t0) is the mapping matrix. The uncertainty of χ0 at the initial moment is P χ , then at the final moment ξ f The uncertainty is P ξ =Φ χ2ξ (t f ,t0)P χ [Φ χ2ξ (t f ,t0)] T (17) Therefore, the present invention directly uses the numerical difference method to calculate Φ χ2ξ (t f ,t0), generate six sets of orthogonal δχ0, corresponding to the right ascension / declination change, right ascension and declination change rate change, distance change, and distance change rate change. After six orbit predictions, six sets of δξ are obtained. f , Φ χ2ξ (t f ,t0) is calculated according to the following formula Therefore, the track point cloud is predicted seven times in parallel. The first time is to obtain the state prediction value, and the remaining six times are to obtain Φ χ2ξ (t f ,t0), and then the predicted orbit uncertainty is obtained.

3. The method for rapid association of optical arc segments of space targets based on point cloud intersection according to claim 1, characterized in that: The specific implementation method of step 4 fast point cloud matching is: In order to achieve fast point cloud intersection calculation, the present invention uses the minimum bounding box method to quickly lock the overlapping area of ​​the two point clouds. The two point clouds formed by the two arc segments are denoted as C1 and C2. The physical quantity X of each data point in C1 is extracted, such as the semi-path p, the improved vernal equinox orbit elements f, g, h, k, and the maximum and minimum values ​​are X respectively. 1max 、X 1min Similarly, extract the X of each data point in C2 and get its maximum and minimum values ​​X 2max 、X 2min , the maximum and minimum values ​​of the minimum bounding box are min{X 1max ,X 2max }、max{X 1min ,X 2min }, obviously, if min{X 1max ,X 2max }<max{X 1min ,X 2min }, then there is no minimum bounding box and the two point clouds do not intersect. If there is a minimum bounding box, discard the data points outside the minimum bounding box in C1 and C2, and you will get X 1min =X 2min 、X 1max =X 2max . After multiple rounds of iterations to find the minimum bounding box for p, f, g, h, k, the number of data points will be significantly reduced, and even converge to a very small area containing the intersection of the two point clouds. Taking into account the observation error, the dynamic model error, and the sampling interval of the distance / distance change rate, the minimum bounding box can be appropriately enlarged. For example, the minimum and maximum values ​​of the minimum bounding box can be taken as min{X 1max ,X 2max }+ΔX、max{X 1min ,X 2min In addition, the present invention does not perform the minimum bounding box operation on the true longitude L, because the true longitude L is a fast variable and is not as stable as other orbital elements. If two arcs belong to the same spatial target, that is, there is a pair of Gaussian samples that "intersect". Whether a pair of Gaussian samples "intersect" can be statistically judged using the Mahalanobis distance, which is defined as Among them, X1 and P1 are the mean and covariance of the first Gaussian sample, and X2 and P2 are the mean and covariance of the second Gaussian sample. If two Gaussian samples intersect, χ 2 It should satisfy the chi-square distribution with 6 degrees of freedom. If there is a pair of Gaussian samples with a Mahalanobis distance less than It is considered that the two arcs may belong to the same spatial target. After obtaining the overlapping area of ​​the two trajectory point clouds using the minimum bounding box operation, each point cloud now contains only a small number of sample points. The remaining sample points are then combined pairwise and the Mahalanobis distance of all combinations is calculated. If a pair of combinations has a Mahalanobis distance less than 4.1, the two arcs are considered to likely belong to the same spatial object.

4. The method for rapid association of optical arc segments of space targets based on point cloud intersection according to claim 1, characterized in that: The specific implementation method of step 5 is to determine the initial trajectory based on point cloud intersection: If the two arc segments obtained after the point cloud intersection are related, the two point clouds after the minimum bounding box operation are written as Gaussian mixture models in, and is the mean and covariance of the i-th Gaussian sample in the first point cloud, and is the mean and covariance of the i-th Gaussian sample in the second point cloud. Taking the first point cloud as prior information and the second point cloud as observation, according to the Bayesian estimation principle, we get in, Further in Finally, Equation (22) can be written as Among them, w ij is the weight coefficient: The initial orbit determination results of the two arc segments are

Images