Method and system for associating space debris observation data with cataloguing orbit data

Abstract

The invention provides a correlation method and system for space debris observation data and cataloguing orbit data. The method comprises the following steps: firstly, within a time range of currently processed batch observation data, recurring to an orbit matrix distributed at equal time intervals through an orbit of an orbit recursion cataloguing target, and generating a position interpolation function based on a time sequence and the orbit matrix; and for observation data in a processing time range, firstly unifying the observation data and an observation station position to a J2000.0 inertial coordinate system, and obtaining a multi-target position prediction value corresponding to observation time through an interpolation function. Calculating the difference between the predicted value and the measured value of the observation data; and calculating an associated object and confidence according to the difference. The method is easy to implement, high in efficiency during batch observation data processing and suitable for association of observation data of different observation devices and different types of orbit space debris and catalogued orbit data.

Description

Technical Field

[0001] The present invention relates to the technical field of space situation awareness, and particularly to a method and system for associating space debris observation data with cataloged orbit data. Background Art

[0002] With the increasing frequency of space activities in recent years, the number of space debris in the Earth's orbit has increased sharply. These space debris, including abandoned satellites, rocket remnants, and debris, pose a serious threat to the spacecraft operating normally in orbit and increase the collision risk between spacecraft. Therefore, it is particularly important to monitor and catalog the orbits of space debris.

[0003] Currently, after obtaining the observation data of space debris through observation equipment, the primary task is to associate these observation data with the cataloged targets. However, observation equipment usually can only obtain limited orbit measurement data and cannot directly obtain the morphology of space debris for identity confirmation. Therefore, associating the observation data with the cataloged orbit data has become a commonly used association method.

[0004] Regarding space debris observation data and cataloged orbital data, the patent (Clustering-based Observation Data and Target Library Element Number Association Method and Apparatus, Patent No. CN 115017386B) discloses a clustering-based observation data and target library element number association method and apparatus. First, based on the space target's observation data, the orbital elements of the space target at a first moment are obtained. Second, the orbital elements of the space target at the first moment are converted into the orbital elements of the space target at a second moment. Then, three parameters are selected from the orbital elements to form two parameter spaces. Finally, two clustering operations are performed in the two parameter spaces on the orbital elements of the space target at the second moment and the orbital elements of the cataloged targets in the preset target library at the second moment, thereby finding the cataloged targets associated with the space target in the preset target library. The patent (A Method for Observation Data-Cataloged Target Association Matching, Patent No. CN 112945182B) discloses an observation data-cataloged target association matching method. A three-layer decision model for the association of observation data with cataloged targets is established. The associated targets are preliminarily screened based on the spatial pointing parameters at the reference time. For the preliminarily screened targets, the observation residuals are calculated separately, and the associated targets are finely screened based on the residual statistical parameters. For the finely screened targets, the association results are selected based on the optimization principle. To address the radar multi-target tracking problem, the literature (Chen Xiang, Research on Data Association Algorithm in the Background of Space Debris, Electronic Measurement Technology, 2017, Vol. 40, No. 8) proposed an improved multi-target data association algorithm suitable for low-orbit space background based on the joint probabilistic data association algorithm. The literature (Huang Qiushi, Short-arc Association Algorithm for Space Targets Based on Sine Fitting, China Space Science and Technology, 2020, Vol. 40, No. 5) associates the space targets observed intermittently by the star sensor, studies the variation law of the right ascension and declination of the new and old tracks, and proposes a short-arc association algorithm for space targets based on sine fitting. Reference (Xu Zhanwei, Probabilistic Data Association Method for Space Target Tracking, Acta Astronomical Sciences, 2017, Vol. 58, No. 3) organically integrates Kalman filtering and probabilistic data association in the optical observation of space targets to form a complete automatic tracking processing chain, realizing a space target adaptive tracking method that integrates multiple measurements within the wave gate. The requirements for the wave gate size can be greatly relaxed, and selecting a larger wave gate can greatly improve the applicability of tracking, especially when the forecast is not accurate enough, it can improve the tracking efficiency. Reference (Lei Xiangxu, Research on Initial Orbit Correlation of Space Debris with Very Short Angle Arc, Acta Geodaetica et Cartographica Sinica, 2021, Vol. 50, No. 2) aims to determine whether two sets of initial orbit parameters come from the same target without any prior information, that is, the initial orbit correlation problem. Based on the analysis of the relationship between orbit parameter error and position error, a purely geometric initial orbit correlation algorithm that does not require initial orbit error information is proposed.

[0005] Therefore, in combination with the current correlation requirements of a large amount of space debris observation data, there is an urgent need to propose a general, fast, and accurate method for correlating observation data and cataloged orbit data applicable to different types of observation equipment and space debris of different orbit types. Summary of the Invention

[0006] Aiming at the defects in the prior art, the purpose of the present invention is to provide a method and system for correlating space debris observation data and cataloged orbit data.

[0007] According to a method for correlating space debris observation data and cataloged orbit data provided by the present invention, the method includes the following steps:

[0008] Step S1: Determine the starting point T s and the ending point T e ;

[0009] Step S2: Recursively calculate the orbits of cataloged targets to the corresponding time series The orbital positions and velocities in the J2000.0 inertial coordinate system, and the result is recorded as the orbital matrix O;

[0010] Step S3: Based on the time series and the orbital matrix O, use the numpy / scipy package to generate a position interpolation function f;

[0011] Step S4: Input the current observation data to be correlated;

[0012] Step S5: Obtain the observation time series The measured values of the observation data in the J2000.0 inertial coordinate system and the station position matrix P;

[0013] Step S6: Generate the position matrix Q of the cataloged target corresponding to the observation time series using the interpolation function f and the observation time series ;

[0014] Step S7: Calculate the predicted values of the observables of the position matrix Q relative to the station position matrix P;

[0015] Step S8: Calculate the difference between the predicted values and the measured values of the observation data;

[0016] Step S9: Calculate the correlation object and confidence level according to the difference;

[0017] Step S10: For the current batch of observation data to be processed, repeat steps S4 to S9 until the end.

[0018] Preferably, in step S2:

[0019] The time series is T s to T e a time series evenly distributed at equal time intervals within the range, which is a one-dimensional vector:

[0020]

[0021] the time series corresponding to the cataloging target the orbit matrix O of is a three-dimensional matrix with dimensions M×6×N, where the 0th dimension corresponds to the time series the length M of the 0th dimension is the length of the time series the length of the 1st dimension 6 corresponds to position and velocity, and the length of the 2nd dimension N corresponds to the number of cataloging targets;

[0022] In the step S5: If the length of the observation time series is K, then the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

[0023] Preferably, the calculation of the predicted value of the observed quantity in the step S7 includes calculating the relative position matrix, the distance matrix, and the relative pointing matrix;

[0024] The calculation of the predicted value of the observed quantity includes the following steps:

[0025] Step S7.1: Calculate the relative position matrix R:

[0026] R = Q - P (Formula 2)

[0027] where, the size of the target position matrix Q is K×3×N, and the size of the station position matrix P is K×3. When subtracting the two, the broadcasting function of numpy arrays needs to be used. During the calculation, P is expanded to K×3×1 for operation, and the size of the resulting relative position matrix R is K×3×N;

[0028] Step S7.2: Calculate the distance matrix ρ. The distance matrix ρ is the norm calculated along the 1st dimension of the relative position matrix R; the size of the distance matrix ρ is K×N;

[0029] Step S7.3: Calculate the relative pointing matrix u:

[0030] u = R / ρ (Formula 3)

[0031] where, the size of the relative position matrix R is K×3×N, and the size of the distance matrix ρ is K×N. When dividing the two, the broadcasting function of numpy arrays needs to be used. During the calculation, ρ is expanded to K×1×N for operation, and the size of the resulting relative pointing matrix u is K×3×N.

[0032] Preferably, the calculation of the difference between the predicted value and the measured value of the ranging type observation quantity in step S8 includes calculating a distance difference matrix and the average distance difference within the observation duration;

[0033] The calculation of the difference between the predicted value and the measured value of the angle measurement type observation quantity includes calculating an angle measurement inner product matrix, an angle measurement difference matrix, and the average angle measurement difference within the observation duration;

[0034] The calculation of the difference between the predicted value and the measured value of the ranging type observation quantity includes the following steps:

[0035] Step S8.1: Calculate the distance difference matrix ρ d :

[0036]

[0037] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. In the calculation, is expanded to K×1 for calculation, and the distance difference matrix ρ d is a two-dimensional matrix with a size of K×N;

[0038] Step S8.2: Calculate the average distance difference within the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ρ d and is a one-dimensional vector with a length of N;

[0039] The calculation of the difference between the predicted value and the measured value of the angle measurement type observation quantity includes the following steps:

[0040] Step S8.3: Calculate the angle measurement inner product matrix u d :

[0041] u d = sum(u·u m ) (Formula 5)

[0042] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. In the calculation, u m is expanded to K×3×1 for element-wise multiplication calculation. The result of u·u m is still a three-dimensional matrix with a size of K×3×N; The sum function represents summing along the 1st dimension of the three-dimensional matrix. The angle measurement inner product matrix u d is a two-dimensional matrix with a size of K×N;

[0043] Step S8.4: Calculate the angle measurement difference matrix a d :

[0044] a d= arccos(u d ) (Formula 6)

[0045] where the arccos function represents the inverse cosine calculation, and the angular difference matrix a d has the same dimensions as the angular inner product matrix u d and is K×N;

[0046] Step S8.5: Calculate the average angular difference within the observation duration which is the result of calculating the mean of the 0th dimension of the distance difference matrix a d and is a one-dimensional vector with a length of N.

[0047] Preferably, the method for calculating associated objects in step S9 includes calculating a difference evaluation factor, splitting the difference evaluation factor to obtain the two smallest values, and judging the associated objects and confidence levels based on these two values;

[0048] If the observed data only has angular data, the evaluation factor is equal to the average angular difference; if the observed data includes ranging and angular values, the evaluation factor is calculated by adjusting the ranging difference and the angular difference through a scaling factor;

[0049] The method for calculating the association confidence level is as follows: if the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the association confidence level is equal to 1; if both the minimum value and the second smallest value are less than the threshold, the association confidence level is calculated according to the formula;

[0050] The method for calculating associated objects includes the following steps:

[0051] Step S9.1: Calculate the difference evaluation factor which is a one-dimensional vector with a length of N;

[0052] Step S9.2: Perform a partition on the difference evaluation factor to obtain the two smallest values;

[0053] Step S9.3: Obtain the minimum value j min and the second smallest value j 2nd ;

[0054] Step S9.4: If the minimum value is greater than the association threshold, the observed data has no matching object; if the minimum value is less than the threshold, the observed data is associated with the cataloged orbit corresponding to the minimum value;

[0055] If the observed data only has angular data, the evaluation factor is equal to the average angular difference If the observed data includes ranging and angular values, the evaluation factor can be calculated by the following formula:

[0056]

[0057] Since radar equipment usually has a relatively high ranging accuracy but a relatively poor angle measurement accuracy, the contribution of the angle measurement data to the evaluation factor is adjusted through the proportionality factor α, where α ≥ 0. When α = 0, the evaluation factor is completely provided by the ranging data;

[0058] If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the association confidence c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation method of the association confidence c is:

[0059]

[0060] The present invention also provides an association system for space debris observation data and cataloged orbit data. The system includes the following modules:

[0061] Module M1: Determine the start time T s and the end time T e ;

[0062] Module M2: Recursively calculate the orbits of the cataloged targets to the corresponding time series of the orbital positions and velocities in the J2000.0 inertial coordinate system, and record the result as the orbit matrix O;

[0063] Module M3: Based on the time series and the orbit matrix O, generate the position interpolation function f using the numpy / scipy package;

[0064] Module M4: Input the current observation data to be associated;

[0065] Module M5: Obtain the observed time series of the measured values of the observation data in the J2000.0 inertial coordinate system and the station position matrix P;

[0066] Module M6: Generate the position matrix Q of the cataloged targets corresponding to the observed time series using the interpolation function f and the observed time series ;

[0067] Module M7: Calculate the predicted values of the observables of the position matrix Q relative to the station position matrix P;

[0068] Module M8: Calculate the difference between the predicted values and the measured values of the observation data;

[0069] Module M9: Calculate the associated objects and confidence levels based on the differences;

[0070] Module M10: Repeatedly call Module M4 to Module M9 for the current batch of observation data to be processed until the end.

[0071] Preferably, in the said Module M2:

[0072] The said time series is T s to T e is a time series with equally spaced time intervals within the range, and is a one-dimensional vector:

[0073]

[0074] The orbit matrix O corresponding to the cataloging target is a three-dimensional matrix with dimensions M×6×N. Among them, the 0th dimension corresponds to the time series The length M of the 0th dimension is the length of the time series The length of the 1st dimension 6 corresponds to position and velocity, and the length of the 2nd dimension N corresponds to the number of cataloging targets;

[0075] In the said Module M5: If the length of the observation time series is K, then the station position

[0076] matrix in the J2000.0 inertial coordinate system has dimensions K×3.

[0077] Preferably, the calculation of the predicted value of the observed quantity in the said Module M7 includes calculating the relative position matrix, distance matrix, and relative pointing matrix;

[0078] The calculation of the predicted value of the observed quantity includes the following modules:

[0079] Module M7.1: Calculate the relative position matrix R:

[0080] R = Q - P (Formula 2)

[0081] Among them, the dimension of the target position matrix Q is K×3×N, and the dimension of the station position matrix P is K×3. When subtracting the two, the broadcasting function of numpy arrays needs to be used, and P is expanded to K×3×1 for operation during the calculation. The dimension of the resulting relative position matrix R is K×3×N;

[0082] Module M7.2: Calculate the distance matrix ρ. The distance matrix ρ is the norm norm calculated along the 1st dimension of the relative position matrix R; the dimension of the distance matrix ρ is K×N;

[0083] Module M7.3: Calculate the relative pointing matrix u:

[0084] u = R / ρ (Formula 3)

[0085] Among them, the relative position matrix R has a size of K×3×N, and the distance matrix ρ has a size of K×N. When dividing the two, the broadcasting function of the numpy array needs to be used. In the calculation, ρ is expanded to K×1×N for operation, and the resulting relative pointing matrix u has a size of K×3×N.

[0086] Preferably, the calculation of the difference between the predicted value and the measured value of the ranging-type observable in the module M8 includes calculating the distance difference matrix and the average distance difference within the observation duration;

[0087] The calculation of the difference between the predicted value and the measured value of the angle-measuring type observable includes calculating the angle-measuring inner product matrix, the angle-measuring difference matrix, and the average angle-measuring difference within the observation duration;

[0088] The calculation of the difference between the predicted value and the measured value of the ranging-type observable includes the following modules:

[0089] Module M8.1: Calculate the distance difference matrix ρ d :

[0090]

[0091] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. In the calculation, is expanded to K×1 for calculation, and the distance difference matrix ρ d is a two-dimensional matrix with a size of K×N;

[0092] Module M8.2: Calculate the average distance difference within the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ρ d and is a one-dimensional vector with a length of N;

[0093] The calculation of the difference between the predicted value and the measured value of the angle-measuring type observable includes the following modules:

[0094] Module M8.3: Calculate the angle-measuring inner product matrix u d :

[0095] u d = sum(u·u m (Formula 5)

[0096] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. In the calculation, u m is expanded to K×3×1 for element-wise multiplication calculation, and the result of u·u m is still a three-dimensional matrix with a size of K×3×N; the sum function represents summing along the 1st dimension of the three-dimensional matrix, and the angle-measuring inner product matrix ud is a two-dimensional matrix with dimensions K×N;

[0097] Module M8.4: Calculate the angular difference matrix a d :

[0098] a d = arccos(u d )(Formula 6)

[0099] where the arccos function represents the inverse cosine calculation, and the angular difference matrix a d has the same dimensions as the angular inner product matrix u d and is K×N;

[0100] Module M8.5: Calculate the average angular difference within the observation duration is the result of calculating the mean of the 0th dimension of the distance difference matrix a d and is a one-dimensional vector with length N.

[0101] Preferably, the associated object calculation system in the module M9 includes calculating a difference evaluation factor, partitioning the difference evaluation factor to obtain the two smallest values, and judging the associated object and confidence level according to these two values;

[0102] If the observed data only has angular measurement data, the evaluation factor is equal to the average angular difference; if the observed data includes ranging and angular measurement values, the evaluation factor is calculated by adjusting the ranging difference and the angular difference through a scaling factor;

[0103] The calculation system for the association confidence level is: if the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the association confidence level is equal to 1; if both the minimum value and the second smallest value are less than the threshold, the association confidence level is calculated according to the formula;

[0104] The associated object calculation system includes the following modules:

[0105] Module M9.1: Calculate the difference evaluation factor is a one-dimensional vector with length N;

[0106] Module M9.2: Partition the difference evaluation factor to obtain the two smallest values;

[0107] Module M9.3: Obtain the minimum value j min and the second smallest value j 2nd ;

[0108] ​

[0109] If the observed data only contains angular measurement data, the evaluation factor is equal to the average angular measurement difference If the observed data includes range measurement and angular measurement values, the evaluation factor can be calculated by the following formula:

[0110]

[0111] Since radar equipment usually has a higher range measurement accuracy and a relatively lower angular measurement accuracy, the contribution of angular measurement data to the evaluation factor is adjusted through the proportionality factor α, where α ≥ 0. When α = 0, that is, the evaluation factor is completely provided by range measurement data;

[0112] If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the association confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation system of the association confidence level c is:

[0113]

[0114] Compared with the prior art, the present invention has the following beneficial effects:

[0115] 1. The present invention can efficiently process a large amount of observed data. By generating orbit recurrence and position interpolation functions, repeated calculations are reduced, and the speed of data association is significantly improved; it is applicable to space debris of different observation equipment and different types of orbits, and is applicable to space debris at different orbit altitudes, showing its wide applicability;

[0116] 2. The present invention calculates the difference between the predicted value and the measured value of the observed data, and calculates the associated object and confidence level based on the difference, ensuring the accuracy of the association result; not only gives the associated object, but also provides the association confidence level, providing strong support for subsequent decision-making and analysis;

[0117] 3. The present invention is based on the Python environment and uses the numpy / scipy package for data processing and calculation. The technical implementation is simple and easy to promote and apply; the method of the present invention is reasonable, the calculation is simple, and the implementation is easy, and it can be effectively applied to the actual association work of space debris observed data and cataloged orbit data. BRIEF DESCRIPTION OF THE DRAWINGS

[0118] By reading the following detailed description of non-limiting embodiments with reference to the accompanying drawings, other features, objects, and advantages of the present invention will become more apparent:

[0119] Figure 1 is the flow chart of the present invention;

[0120] Figure 2 Schematic diagram of the associated time consumption for multiple arc segments;

[0121] Figure 3 Schematic diagram of the associated confidence for multiple arc segments. Detailed implementation manner

[0122] The present invention will be described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several changes and improvements can still be made. These all belong to the protection scope of the present invention.

[0123] Example 1:

[0124] Referring to Figure 1 and Figure 2 , according to a method for associating space debris observation data with cataloged orbit data provided by the present invention, the method includes the following steps:

[0125] Step S1: Determine the starting point T s and the ending point T e ;

[0126] Step S2: Propagate the orbit of the cataloged target to the corresponding time sequence The orbit position and velocity in the J2000.0 inertial coordinate system, and the result is recorded as the orbit matrix O; the time sequence is a time sequence evenly distributed at equal time intervals within the range from T s to T e , and is a one-dimensional vector:

[0127]

[0128] The orbit matrix O of the cataloged target corresponding to the time sequence is a three-dimensional matrix with a size of M×6×N, where the 0th dimension corresponds to the time sequence The length M of the 0th dimension is the length of the time sequence , the length 6 of the 1st dimension corresponds to the position and velocity, and the length N of the 2nd dimension corresponds to the number of cataloged targets;

[0129] Step S3: Based on the time sequence and the orbit matrix O, generate a position interpolation function f using the numpy / scipy package;

[0130] Step S4: Input the current observation data to be associated;

[0131] Step S5: Obtain the observed time series The measured values of the observed data in the J2000.0 inertial coordinate system and the station position matrix P; if the length of the observed time series is K, then the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

[0132] Step S6: Generate the position matrix Q of the cataloging target corresponding to the observed time series from the interpolation function f and the observed time series ;

[0133] Step S7: Calculate the predicted values of the observables of the position matrix Q relative to the station position matrix P; the calculation of the predicted values of the observables includes calculating the relative position matrix, the distance matrix, and the relative pointing matrix;

[0134] The calculation of the predicted values of the observables includes the following steps:

[0135] Step S7.1: Calculate the relative position matrix R:

[0136] R = Q - P (Formula 2)

[0137] where the size of the target position matrix Q is K×3×N, and the size of the station position matrix P is K×3. When subtracting the two, the broadcasting function of the numpy array needs to be used, and P is expanded to K×3×1 for calculation during the calculation. The size of the resulting relative position matrix R is K×3×N;

[0138] Step S7.2: Calculate the distance matrix ρ. The distance matrix ρ is the norm calculated along the first dimension of the relative position matrix R; the size of the distance matrix ρ is K×N;

[0139] Step S7.3: Calculate the relative pointing matrix u:

[0140] u = R / ρ (Formula 3)

[0141] where the size of the relative position matrix R is K×3×N, and the size of the distance matrix ρ is K×N. When dividing the two, the broadcasting function of the numpy array needs to be used, and ρ is expanded to K×1×N for calculation during the calculation. The size of the resulting relative pointing matrix u is K×3×N.

[0142] Step S8: Calculate the difference between the predicted values and the measured values of the observed data; the calculation of the difference between the predicted values and the measured values of the ranging observables includes calculating the distance difference matrix and the average distance difference within the observation duration;

[0143] The calculation of the difference between the predicted values and the measured values of the angle-measuring observables includes calculating the angle inner product matrix, the angle difference matrix, and the average angle difference within the observation duration;

[0144] The calculation of the difference between the predicted value and the measured value of the ranging-type observation quantity includes the following steps:

[0145] Step S8.1: Calculate the distance difference matrix ρd:

[0146]

[0147] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. During the calculation, is extended to K×1 for calculation. The calculated distance difference matrix ρd is a two-dimensional matrix with a size of K×N;

[0148] Step S8.2: Calculate the average distance difference during the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ρ d and is a one-dimensional vector with a length of N;

[0149] The calculation of the difference between the predicted value and the measured value of the angle measurement-type observation quantity includes the following steps:

[0150] Step S8.3: Calculate the angle measurement inner product matrix u d :

[0151] u d = sum(u·u m )(Formula 5)

[0152] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. During the calculation, u m is extended to K×3×1 for element-by-element multiplication calculation. The result of u·u m is still a three-dimensional matrix with a size of K×3×N; the sum function represents summing along the 1st dimension of the three-dimensional matrix. The angle measurement inner product matrix u d is a two-dimensional matrix with a size of K×N;

[0153] Step S8.4: Calculate the angle measurement difference matrix a d :

[0154] a d = arccos(u d )(Formula 6)

[0155] Among them, the arccos function represents the inverse cosine calculation. The size of the angle measurement difference matrix a d is the same as that of the angle measurement inner product matrix u d and is K×N;

[0156] Step S8.5: Calculate the average angle measurement difference during the observation duration The result of calculating the mean value of the 0th dimension of the distance difference matrix a d is a one-dimensional vector with a length of N.

[0157] Step S9: Calculate the associated object and confidence according to the difference; the method for calculating the associated object includes calculating the difference evaluation factor, splitting the difference evaluation factor to obtain the two smallest values, and judging the associated object and confidence according to these two values;

[0158] If the observed data only has angle measurement data, the evaluation factor is equal to the average angle measurement difference; if the observed data includes distance measurement and angle measurement values, the evaluation factor is calculated by adjusting the distance measurement difference and the angle measurement difference through a proportionality factor;

[0159] Refer to Figure 3 , the method for calculating the association confidence is: if the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the association confidence is equal to 1; if both the minimum value and the second smallest value are less than the threshold, calculate the association confidence according to the formula;

[0160] The method for calculating the associated object includes the following steps:

[0161] Step S9.1: Calculate the difference evaluation factor is a one-dimensional vector with a length of N;

[0162] Step S9.2: Perform a partition on the difference evaluation factor to obtain the two smallest values;

[0163] Step S9.3: Obtain the minimum value j min and the second smallest value j 2nd ;

[0164] Step S9.4: If the minimum value is greater than the association threshold, there is no matching object for this observed data; if the minimum value is less than the threshold, this observed data is associated with the cataloged orbit corresponding to the minimum value;

[0165] If the observed data only has angle measurement data, the evaluation factor is equal to the average angle measurement difference If the observed data includes distance measurement and angle measurement values, the evaluation factor can be calculated by the following formula:

[0166]

[0167] Since radar equipment usually has a higher ranging accuracy and a relatively lower angle measurement accuracy, the contribution of the angle measurement data to the evaluation factor is adjusted through a proportionality factor α, α≥0. When α = 0, that is, the evaluation factor is completely provided by the ranging data;

[0168] If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the associated confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation method of the associated confidence level c is as follows:

[0169]

[0170] Step S10: For the current batch of observation data to be processed, repeat Step S4 to Step S9 until the end.

[0171] The present invention also provides a correlation system for space debris observation data and cataloged orbit data. The correlation system for space debris observation data and cataloged orbit data can be implemented by executing the process steps of the correlation method for space debris observation data and cataloged orbit data. That is, those skilled in the art can understand the correlation method for space debris observation data and cataloged orbit data as a preferred implementation manner of the correlation system for space debris observation data and cataloged orbit data.

[0172] Example 2:

[0173] The present invention also provides a correlation system for space debris observation data and cataloged orbit data. The system includes the following modules:

[0174] Module M1: Determine the start time T s and the end time T e ;

[0175] Module M2: Recursively calculate the orbits of the cataloged targets to the orbital positions and velocities at the corresponding moment sequences in the J2000.0 inertial coordinate system. The result is recorded as the orbit matrix O; the moment sequence is s from T e to T

[0176]

[0177] The orbit matrix O of the cataloged target corresponding to the moment sequence is a three-dimensional matrix with a size of M×6×N. Among them, the 0th dimension corresponds to the moment sequence The length M of the 0th dimension is the length of the moment sequence The length of the 1st dimension is 6 corresponding to the position and velocity, and the length of the 2nd dimension is N corresponding to the number of cataloged targets;

[0178] Module M3: Based on the moment sequence Using the orbital matrix O, generate the position interpolation function f with the numpy / scipy package;

[0179] Module M4: Input the current observation data to be associated;

[0180] Module M5: Obtain the observation time series The measured values of the observation data in the J2000.0 inertial coordinate system and the station position matrix P; if the length of the observation time series is K, then the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

[0181] Module M6: Generate the position matrix Q of the catalog target corresponding to the observation time series from the interpolation function f and the observation time series ;

[0182] Module M7: Calculate the predicted value of the observable relative to the station position matrix P; the calculation of the predicted value of the observable includes calculating the relative position matrix, the distance matrix, and the relative pointing matrix;

[0183] The calculation of the predicted value of the observable includes the following modules:

[0184] Module M7.1: Calculate the relative position matrix R:

[0185] R = Q - P (Formula 2)

[0186] where, the size of the target position matrix Q is K×3×N, the size of the station position matrix P is K×3, and the subtraction of the two requires the use of the broadcasting function of the numpy array. In the calculation, P is expanded to K×3×1 for operation, and the size of the resulting relative position matrix R is K×3×N;

[0187] Module M7.2: Calculate the distance matrix ρ. The distance matrix ρ is the norm calculated along the first dimension of the relative position matrix R; the size of the distance matrix ρ is K×N;

[0188] Module M7.3: Calculate the relative pointing matrix u:

[0189] u = R / ρ (Formula 3)

[0190] where, the size of the relative position matrix R is K×3×N, the size of the distance matrix ρ is K×N, and the division of the two requires the use of the broadcasting function of the numpy array. In the calculation, ρ is expanded to K×1×N for operation, and the size of the resulting relative pointing matrix u is K×3×N.

[0191] Module M8: Calculate the difference between the predicted value and the measured value of the observation data; the calculation of the difference between the predicted value and the measured value of the ranging observable includes calculating the distance difference matrix and the average distance difference within the observation duration;

[0192] The calculation of the difference between the predicted value and the measured value of the angle measurement type observation quantity includes calculating the angle measurement inner product matrix, the angle measurement difference matrix, and the average angle measurement difference within the observation duration;

[0193] The calculation of the difference between the predicted value and the measured value of the distance measurement type observation quantity includes the following modules:

[0194] Module M8.1: Calculate the distance difference matrix ρd:

[0195]

[0196] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. During the calculation, is extended to K×1 for calculation, and the distance difference matrix ρ is calculated. d is a two-dimensional matrix with a size of K×N;

[0197] Module M8.2: Calculate the average distance difference within the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ρ d and is a one-dimensional vector with a length of N;

[0198] The calculation of the difference between the predicted value and the measured value of the angle measurement type observation quantity includes the following modules:

[0199] Module M8.3: Calculate the angle measurement inner product matrix u d :

[0200] u d = sum(u·u m ) (Formula 5)

[0201] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. During the calculation, u m is extended to K×3×1 for element-wise multiplication calculation, and the result of u·u m is still a three-dimensional matrix with a size of K×3×N; the sum function represents summing along the 1st dimension of the three-dimensional matrix, and the angle measurement inner product matrix u d is a two-dimensional matrix with a size of K×N;

[0202] Module M8.4: Calculate the angle measurement difference matrix ad:

[0203] a d = arccos(u d ) (Formula 6)

[0204] Among them, the arccos function represents the inverse cosine calculation, and the angle measurement difference matrix ad The size is the same as that of the angular measurement inner product matrix u d and is K×N;

[0205] Module M8.5: Calculate the average angular measurement difference within the observation duration It is the result of calculating the mean value of the 0th dimension of the distance difference matrix a d and is a one-dimensional vector with a length of N.

[0206] Module M9: Calculate the associated object and confidence level based on the difference; The associated object calculation system includes calculating the difference evaluation factor, partitioning the difference evaluation factor to obtain the two smallest values, and judging the associated object and confidence level based on these two values;

[0207] If the observation data only contains angular measurement data, the evaluation factor is equal to the average angular measurement difference; If the observation data includes distance measurement and angular measurement values, the evaluation factor is calculated by adjusting the distance measurement difference and the angular measurement difference through a proportionality factor;

[0208] The calculation system for the associated confidence level is: If the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the associated confidence level is equal to 1; If both the minimum value and the second smallest value are less than the threshold, the associated confidence level is calculated according to the formula;

[0209] The associated object calculation system includes the following modules:

[0210] Module M9.1: Calculate the difference evaluation factor It is a one-dimensional vector with a length of N;

[0211] Module M9.2: Partition the difference evaluation factor to obtain the two smallest values;

[0212] Module M9.3: Obtain the minimum value j min and the second smallest value j 2nd ;

[0213] Module M9.4: If the minimum value is greater than the associated threshold, there is no matching object for this observation data; If the minimum value is less than the threshold, this observation data is associated with the cataloged orbit corresponding to the minimum value;

[0214] If the observation data only contains angular measurement data, the evaluation factor is equal to the average angular measurement difference If the observation data includes distance measurement and angular measurement values, the evaluation factor can be calculated by the following formula:

[0215]

[0216] Since radar equipment usually has high ranging accuracy but relatively poor angle measurement accuracy, the contribution of angle measurement data to the evaluation factor is adjusted through a scale factor α, where α ≥ 0. When α = 0, the evaluation factor is completely provided by the ranging data;

[0217] If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the association confidence c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation system for the association confidence c is:

[0218]

[0219] Module M10: Repeatedly call Module M4 to Module M9 for the current batch of observation data to be processed until the end.

[0220] Example 3:

[0221] A fast association method for space debris observation data and cataloged orbit data includes the following steps:

[0222] Step 1, determine the start time T s and the end time T e of the time range of the current batch of observation data to be processed, that is, the current processed observation data are all within the range not exceeding T s to T e range.

[0223] Step 2, recursively calculate the orbits of the cataloged targets to the corresponding time series of the orbital positions and velocities in the J2000.0 inertial coordinate system. The result is recorded as the orbit matrix O, where the time series is a time series evenly distributed at equal time intervals within the range of T s to T e and is a one-dimensional vector:

[0224]

[0225] The orbit matrix O of the cataloged target corresponding to the time series is a three-dimensional matrix with dimensions M×6×N, where the 0th dimension corresponds to the time series The length M of the 0th dimension is the length of the time series The length of the 1st dimension is 6 corresponding to the position and velocity, and the length of the 2nd dimension is N corresponding to the number of cataloged targets.

[0226] Step 3, based on the time series and the orbit matrix O, generate a position interpolation function f using the numpy / scipy package.

[0227] Step 4: Input the current observation data to be associated. Typical observation data includes the angle-only observation data of optical devices and the range and angle observation data of radar devices.

[0228] Step 5: Obtain the observation time series The measured values of the observation data in the J2000.0 inertial coordinate system and the station position matrix P. If the length of the observation time series is K, then the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

[0229] Step 6: Generate the position matrix Q of the catalog target corresponding to the observation time series from the interpolation function f and the observation time series

[0230] The position matrix Q of the catalog target corresponding to the observation time series is a three-dimensional matrix with a size of K×3×N. The 0th dimension corresponds to the observation time series The length K of the 0th dimension is the length of the observation time series The length 3 of the 1st dimension corresponds to the position, and the length N of the 2nd dimension corresponds to the number of catalog targets.

[0231] Step 7: Calculate the predicted value of the observable quantity of the position matrix Q relative to the station position matrix P.

[0232] Step 8: Calculate the difference between the predicted value and the measured value of the observation data.

[0233] Step 9: Calculate the associated object and confidence according to the difference.

[0234] Step 10: For the batch of observation data to be processed currently, repeat Steps 4 to 9 until the end.

[0235] The measured values of the observation data in the J2000.0 inertial coordinate system in Step 5 are divided into range measurement type data and angle measurement type data according to different observation devices. Among them, typical optical devices only obtain angle measurement data, and typical radar devices can obtain range and angle measurement data. The range data sequence corresponding to the observation time series is denoted as and is a one-dimensional vector with a length of K. The angle measurement data corresponding to the observation time series is converted into the pointing matrix u m in the J2000.0 inertial coordinate system, with a size of K×3. The 0th dimension of the pointing matrix u corresponds to the observation time series The length K of the 0th dimension is the length of the observation time series The pointing vector at any moment (with a size of 3, corresponding to the pointing matrix u m ​​The first dimension) are all unit vectors, representing the measured values of the pointing vector of the space debris relative to the measurement station in the J2000.0 inertial coordinate system.

[0236] In step 7, the calculation of the predicted value of the observable contains the following steps:

[0237] Step 7.1, calculate the relative position matrix R.

[0238] R = Q - P (Formula 2)

[0239] Among them, the target position matrix Q (size K×3×N), the measurement station position matrix P (size K×3), and the subtraction of the two requires the use of the broadcasting function of the numpy array. During the calculation, P is expanded to K×3×1 for operation. The size of the resulting relative position matrix R is K×3×N.

[0240] Step 7.2, calculate the distance matrix ρ. The distance matrix ρ is the norm calculated along the first dimension of the relative position matrix R. The size of the distance matrix ρ is K×N.

[0241] Step 7.3, calculate the relative pointing matrix u.

[0242] u = R / ρ (Formula 3)

[0243] Among them, the relative position matrix R has a size of K×3×N, and the distance matrix ρ has a size of K×N. The division of the two requires the use of the broadcasting function of the numpy array. During the calculation, ρ is expanded to K×1×N for operation. The size of the resulting relative pointing matrix u is K×3×N.

[0244] In step 8, the calculation of the difference between the predicted value and the measured value of the ranging observable contains the following steps:

[0245] Step 8.1, calculate the distance difference matrix ρd.

[0246]

[0247] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. During the calculation, is expanded to K×1 for calculation. Calculate the distance difference matrix ρ d is a two-dimensional matrix with a size of K×N.

[0248] Step 8.2, calculate the average distance difference during the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ρ d and is a one-dimensional vector with a length of N.

[0249] In step 8, the calculation of the difference between the predicted value and the measured value of the angle measurement type observation quantity includes the following steps:

[0250] Step 8.3, calculate the angle measurement inner product matrix u d .

[0251] u d = sum(u·u m ) (Formula 5)

[0252] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. During the calculation, u m is expanded to K×3×1 for element-by-element multiplication calculation. The result of u·u m is still a three-dimensional matrix with a size of K×3×N. The sum function represents summation along the first dimension of the three-dimensional matrix. The angle measurement inner product matrix u d is a two-dimensional matrix with a size of K×N.

[0253] Step 8.4, calculate the angle measurement difference matrix ad.

[0254] a d = arccos(u d ) (Formula 6)

[0255] Among them, the arccos function represents inverse cosine calculation. The angle measurement difference matrix a d has the same size as the angle measurement inner product matrix u d , which is K×N.

[0256] Step 8.5, calculate the average angle measurement difference during the observation duration is the result of calculating the mean value of the 0th dimension of the distance difference matrix ad, which is a one-dimensional vector with a length of N.

[0257] In step 9, the calculation method of the associated object includes the following steps:

[0258] Step 9.1, calculate the difference evaluation factor ( is a one-dimensional vector with a length of N).

[0259] Step 9.2, partition the difference evaluation factor to obtain the two smallest values.

[0260] Step 9.3, obtain the minimum value j min and the second minimum value j 2nd .

[0261] Step 9.4, if the minimum value is greater than the correlation threshold, there is no matching object for this observation data. If the minimum value is less than the threshold, this observation data is associated with the cataloged orbit corresponding to the minimum value.

[0262] In Step 9, if the observation data only has angular measurement data, the evaluation factor is equal to the average angular measurement difference If the observation data includes range measurement and angular measurement values, the evaluation factor can be calculated by the following formula

[0263]

[0264] Since radar equipment usually has a high range measurement accuracy, while the angular measurement accuracy is relatively poor, the contribution of angular measurement data to the evaluation factor can be adjusted by a scaling factor α (α≥0). When α = 0, that is, the evaluation factor is completely provided by range measurement data.

[0265] In Step 9, if the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the association confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation method of the association confidence level c is

[0266]

[0267] The technical problem to be solved by the present invention is to provide a fast association method for space debris observation data and cataloged orbit data, which can be applicable to space debris of different observation equipment and different orbit types, accurately associate and give the association confidence level.

[0268] The present invention provides a method that can quickly implement the association of space debris observation data and cataloged orbit data. Based on the python environment, the association processing flow of batch observation data based on the numpy / scipy package is given. The difference evaluation factors for different types of observation data and the calculation method of the association confidence level are constructed. The method of the present invention is reasonable, simple in calculation, and easy to implement, and can be effectively applied to space debris at different orbit altitudes, and is applicable to the fast association of observation data of optical and radar equipment with cataloged orbits.

[0269] After obtaining the orbit measurement data of space debris, it is first necessary to associate and match it with the cataloged targets to confirm whether it is a certain cataloged target. With the gradual improvement of current observation equipment, especially space-based observation equipment, which can effectively improve the observation coverage, the orbit measurement data that needs to be processed currently has increased sharply, and it is necessary to complete the matching quickly within a short time. The typical observation data of space debris is divided into angle-only observation data of optical equipment and range and angle observation data of radar equipment. When associating with the orbit information to be cataloged and extrapolating the cataloged orbit, in order to save processing time, first extrapolate the orbits of all cataloged targets to a certain time range, which should cover the batch of observation data currently being processed. If the start time T s , end time T e , then all the observation data currently being processed does not exceed the range from T s to T e .

[0270] Extrapolate the orbits of the cataloged targets to the corresponding time sequence The orbital positions and velocities in the J2000.0 inertial coordinate system, and the result is recorded as the orbit matrix O, where the time sequence is the time sequence evenly distributed at equal time intervals within the range from T s to T e , and it is a one-dimensional vector:

[0271]

[0272] The orbit matrix O corresponding to the time sequence of the cataloged targets is a three-dimensional matrix with a size of M×6×N, where the 0th dimension corresponds to the time sequence The length M of the 0th dimension is the length of the time sequence , the length 6 of the 1st dimension corresponds to the position and velocity, and the length N of the 2nd dimension corresponds to the number of cataloged targets. In the python environment, O is a three-dimensional array of numpy.

[0273] Since the orbit matrix O corresponds to an equally spaced time sequence , and in the association, it is necessary to obtain the positions of multiple targets at the corresponding observation times. Therefore, based on the time sequence and the orbit matrix O, use the numpy / scipy package to generate the position interpolation function f. After obtaining the interpolation function f, for the observation data within the subsequent time range from T s to T e , the interpolation function does not need to be obtained again.

[0274] Input the current observation data to be associated. The data formats of each station equipment may vary, and here they can be uniformly converted to the J2000.0 inertial coordinate system for processing. Obtain the observation time sequence The measured values of the observation data in the J2000.0 inertial coordinate system and the station position matrix P. If the length of the observation time series is K, the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

[0275] The measured values of the observation data in the J2000.0 inertial coordinate system are divided into ranging data and angle measurement data according to different observation devices. Among them, typical optical devices only obtain angle measurement data, and typical radar devices can obtain ranging and angle measurement data. Corresponding to the observation time series The ranging data sequence is denoted as which is a one-dimensional vector with a length of K. Corresponding to the observation time series The angle measurement data is converted into the pointing matrix u m in the J2000.0 inertial coordinate system, with a size of K×3. The 0th dimension of the pointing matrix u corresponds to the observation time series The 0th dimension length K is the length of the observation time series At any moment, the pointing vector (with a size of 3, corresponding to the 1st dimension of the pointing matrix u m ) is a unit vector, representing the measured value of the pointing vector of the space debris relative to the station in the J2000.0 inertial coordinate system.

[0276] The position matrix Q of the cataloged target corresponding to the observation time series is generated by the interpolation function f and the observation time series .

[0277] The position matrix Q of the cataloged target corresponding to the observation time series is a three-dimensional matrix with a size of K×3×N, where the 0th dimension corresponds to the observation time series The 0th dimension length K is the length of the observation time series , the 1st dimension length 3 corresponds to the position, and the 2nd dimension length N corresponds to the number of cataloged targets.

[0278] The 0th dimension of both the position matrix Q and the station position matrix P corresponds to the observation time, so the relative position in the J2000.0 inertial coordinate can be obtained by direct subtraction, that is, the relative position matrix R is

[0279] R = Q - P (Formula 2)

[0280] Among them, the target position matrix Q (size K×3×N), the station position matrix P (size K×3). When subtracting the two, the broadcasting function of the numpy array needs to be used, and P is extended to K×3×1 for operation in the calculation. The size of the resulting relative position matrix R is K×3×N.

[0281] The modulus value of the relative position corresponds to the relative distance, and the normalized vector corresponds to the relative angle measurement. Therefore, the distance matrix ρ is the norm calculated along the first dimension of the relative position matrix R. The size of the distance matrix ρ is K×N.

[0282] Normalize the first dimension of the matrix R using the relative distance matrix ρ to obtain the relative pointing matrix u.

[0283] u = R / ρ (Formula 3)

[0284] Among them, the size of the relative position matrix R is K×3×N, and the size of the distance matrix ρ is K×N. When dividing the two, the broadcasting function of the numpy array needs to be used, and ρ is expanded to K×1×N for calculation during the calculation. The size of the resulting relative pointing matrix u is K×3×N.

[0285] Based on the difference between the predicted value and the measured value of the observed data, the matching object can be judged, and different calculation methods are performed for the ranging-type observed data and the angle-measuring-type observed data.

[0286] The difference between the predicted value and the measured value of the ranging-type observable is first the distance difference matrix ρ d .

[0287]

[0288] Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K, and during the calculation, is expanded to K×1 for calculation. The distance difference matrix ρ d is a two-dimensional matrix with a size of K×N.

[0289] To facilitate the judgment, the distance difference matrix ρ d is averaged along the 0th dimension to obtain the average distance difference within the observation duration is a one-dimensional vector with a length of N.

[0290] The difference between the predicted value and the measured value of the angle-measuring-type observable is first to calculate the angular inner product matrix u d .

[0291] u d = sum(u·u m ) (Formula 5)

[0292] Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. During the calculation, u m is expanded to K×3×1 for element-wise multiplication calculation, u·u mThe result is still a three-dimensional matrix with dimensions K×3×N. The sum function represents summing along the first dimension of the three-dimensional matrix, and the angular inner product matrix u d is a two-dimensional matrix with dimensions K×N.

[0293] The inner product of two unit vectors corresponds to the cosine value of the angle between the two unit vectors. That is, the angular difference matrix ad satisfies

[0294] a d = arccos(u d ) (Formula 6)

[0295] where the arccos function represents the inverse cosine calculation. The angular difference matrix a d has the same dimensions as the angular inner product matrix u d and is K×N.

[0296] Similarly, the average angular difference within the observation duration can be calculated which is the result of calculating the mean value of the 0th dimension of the distance difference matrix a d and is a one-dimensional vector with a length of N.

[0297] Obtain the average distance difference or the average angular difference After that, calculate the difference evaluation factor. If the observed data only contains angular data, the evaluation factor is equal to the average angular difference If the observed data includes ranging and angular values, the evaluation factor can be calculated by the following formula

[0298]

[0299] Since radar equipment usually has a higher ranging accuracy and a relatively lower angular accuracy, the contribution of the angular data to the evaluation factor can be adjusted by the scaling factor α (α≥0). When α = 0, that is, the evaluation factor is completely provided by the ranging data.

[0300] The evaluation factor is a one-dimensional vector with a length of N, corresponding to N cataloged targets. The smaller the evaluation factor value, the smaller the difference between the predicted value and the measured value. Conventional matching considers the minimum value of the evaluation factor as the most likely matching object, and the present invention combines a threshold and optimal matching for association judgment.

[0301] When there is exactly one target within the threshold, its association has a relatively high confidence; when there are multiple targets within the threshold and the difference between the evaluation factors corresponding to the optimal match and the sub-optimal match is not significant, the confidence of associating with the optimal match target is relatively small at this time. When no target is within the threshold, it returns that no cataloged target is associated.

[0302] The calculation of the confidence level requires obtaining the evaluation factors with the minimum value j min and the second minimum value j 2nd . In the Python environment, to obtain the two smallest values, the argpartition interface of the numpy package can be used.

[0303] If the minimum value j min is less than the threshold while the second minimum value j 2nd is greater than the threshold, the associated confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation method of the associated confidence level c is

[0304]

[0305] The effectiveness of the method of the present invention is illustrated below by processing the measured data of space debris. To verify the correctness and speed of the association, the verification data uses arcs with known true values of the associated objects, and only their measured data is used to perform the association, and the association results are compared with the true values. The number of arcs in the verification data is 110. In the ThinkStation P350 Python 3.8 environment, for 20,954 cataloged orbits, it takes 839.7 s to complete steps 1 to 3, and there is no need to repeat the operation for each arc after steps 1 to 3. The time consumption of each observation data in steps 4 to 9 is as shown in the appendix Figure 2 as shown, and the time consumption does not exceed 2 s, with an average time consumption of 0.4 s for each observation data. The association results of the 110 observation data are all correctly associated with the true values, and the output confidence levels are as shown in the appendix Figure 3 as shown, and the association confidence levels of the 110 observation data in the verification data are all greater than 0.944.

[0306] Those skilled in the art can understand this embodiment as a more specific description of Embodiment 1 and Embodiment 2.

[0307] Those skilled in the art know that in addition to implementing the system and its various devices, modules, and units provided by the present invention in the form of pure computer-readable program code, the method steps can be logically programmed to enable the system and its various devices, modules, and units provided by the present invention to be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers, etc. to achieve the same functions. Therefore, the system and its various devices, modules, and units provided by the present invention can be regarded as a kind of hardware component, and the devices, modules, and units included therein for implementing various functions can also be regarded as the structures within the hardware component; the devices, modules, and units for implementing various functions can also be regarded as both software modules for implementing the method and the structures within the hardware component.

[0308] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention. Without conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. A method for associating space debris observation data with cataloged orbit data, characterized in that, The method includes the following steps: Step S1: Determine the start point T of the time range of the batch observation data currently being processed s and the end point T e ; Step S2: Recursively derive the orbit of the cataloging target to the corresponding time series The orbital position and velocity in the J2000.0 inertial coordinate system, and the result is recorded as the orbit matrix O; Step S3: Based on the time series and the orbit matrix O, use the numpy / scipy package to generate a position interpolation function f; Step S4: Input the current observation data to be associated; Step S5: Obtain the observation time series The measured values of the observed data in the J2000.0 inertial coordinate system and the station position matrix P; Step S6: Generate the position matrix Q of the cataloging target corresponding to the observation time series from the interpolation function f and the observed time series of the observed time series​ Step S7: Calculate the predicted value of the observable quantity of the position matrix Q relative to the station position matrix P; Step S8: Calculate the difference between the predicted value and the measured value of the observation data; Step S9: Calculate the associated object and confidence level according to the difference; Step S10: For the batch of observation data to be processed currently, repeat Step S4 to Step S9 until the end.

2. The correlation method between space debris observation data and cataloged orbit data according to claim 1, wherein, In the said Step S2: The said time series is T s to T e The time series evenly distributed at equal time intervals within the range is a one-dimensional vector: Catalog target corresponding time series The orbit matrix O of is a three-dimensional matrix with dimensions M×6×N, where the 0th dimension corresponds to the time series The length M of the 0th dimension is the length of the time series The length of the 1st dimension 6 corresponds to position and velocity, and the length of the 2nd dimension N corresponds to the number of catalog targets; In the step S5: If the length of the observed time series is K, the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

3. The correlation method for space debris observation data and cataloged orbit data according to claim 1, characterized in that In the calculation of the predicted value of the observable quantity in Step S7, it includes calculating the relative position matrix, distance matrix, and relative pointing matrix; The calculation of the predicted value of the observable quantity includes the following steps: Step S7.1: Calculate the relative position matrix R: R = Q - P (Formula 2) Wherein, the size of the target position matrix Q is K×3×N, the size of the station position matrix P is K×3, and when subtracting the two, the broadcast function of the numpy array needs to be used. In the calculation, P is expanded to K×3×1 for operation, and the size of the resulting relative position matrix R is K×3×N; Step S7.2: Calculate the distance matrix ρ. The distance matrix ρ is to calculate the norm along the first dimension of the relative position matrix R; the size of the distance matrix ρ is K×N; Step S7.3: Calculate the relative pointing matrix u: u = R / ρ (Formula 3) Wherein, the size of the relative position matrix R is K×3×N, the size of the distance matrix ρ is K×N, and when dividing the two, the broadcast function of the numpy array needs to be used. In the calculation, ρ is expanded to K×1×N for operation, and the size of the resulting relative pointing matrix u is K×3×N.

4. The correlation method for space debris observation data and cataloged orbit data according to claim 1, characterized in that, In the calculation of the difference between the predicted value and the measured value of the ranging observable quantity in Step S8, it includes calculating the distance difference matrix and the average distance difference within the observation duration; The calculation of the difference between the predicted value and the measured value of the angle measurement observable quantity includes calculating the angle measurement inner product matrix, angle measurement difference matrix, and the average angle measurement difference within the observation duration; The calculation of the difference between the predicted value and the measured value of the ranging observable quantity includes the following steps: Step S8.1: Calculate the distance difference matrix ρ d : Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N, is a one-dimensional vector with a length of K. During the calculation, is expanded to K×1 for calculation, and the distance difference matrix ρ d is a two-dimensional matrix with a size of K×N; Step S8.2: Calculate the average distance difference within the observation duration is the result of calculating the mean of the 0th dimension of the distance difference matrix ρ d which is a one-dimensional vector with a length of N; The calculation of the difference between the predicted value and the measured value of the angle measurement observable quantity includes the following steps: Step S8.3: Calculate the angular inner product matrix u d : u d = sum(u·u m ) (Equation 5) Among them, the relative pointing matrix u is a three-dimensional matrix with a size of K×3×N, which is the measurement pointing matrix u m , with a size of K×3. During the calculation, u m is extended to K×3×1 for element-wise multiplication calculation. The result of u·u m is still a three-dimensional matrix with a size of K×3×N; the sum function represents summation along the first dimension of the three-dimensional matrix. The angular inner product matrix u d is a two-dimensional matrix with a size of K×N; Step S8.4: Calculate the angular difference matrix a d : a d = arccos(u d ) (Formula 6) wherein, the arccos function represents the arccosine calculation, and the angle measurement difference matrix a d has the same size as the angle measurement inner product matrix u d which is K×N; Step S8.5: Calculate the average angular measurement difference within the observation duration is the result of calculating the mean of the 0th dimension of the distance difference matrix a d which is a one-dimensional vector with a length of N.

5. The correlation method for space debris observation data and cataloged orbit data according to claim 1, characterized in that In the said Step S9, the method for calculating the associated object includes calculating the difference evaluation factor, and partitioning the difference evaluation factor to obtain the two smallest values, and judging the associated object and confidence level according to these two values; If the observation data only has angle measurement data, the evaluation factor is equal to the average angle measurement difference; If the observation data includes ranging and angle measurement values, the evaluation factor is calculated by adjusting the ranging difference and angle measurement difference through a proportionality factor; The calculation method of the association confidence level is: if the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the association confidence level is equal to 1; if both the minimum value and the second smallest value are less than the threshold, the association confidence level is calculated according to the formula; The method for calculating the associated object includes the following steps: Step S9.1: Calculate the difference evaluation factor is a one-dimensional vector with a length of N; Step S9.2: Partition the difference evaluation factor to obtain the two smallest values; Step S9.3: Obtain the minimum value j of the difference evaluation factor min and the second minimum value j 2nd ; Step S9.4: If the minimum value is greater than the association threshold, there is no matching object for this observation data; If the minimum value is less than the threshold, this observation data is associated with the cataloged orbit corresponding to the minimum value; If the observed data only contains angular measurement data, the evaluation factor is equal to the average angular measurement difference If the observed data contains range measurement and angular measurement values, the evaluation factor can be calculated by the following formula: Since radar equipment usually has a relatively high ranging accuracy and a relatively poor angle measurement accuracy, the contribution of the angle measurement data to the evaluation factor is adjusted by a scaling factor α, where α ≥ 0. When α = 0, the evaluation factor is completely provided by the ranging data. If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, then the associated confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation method of the associated confidence level c is as follows:

6. A correlation system for space debris observation data and cataloged orbit data, characterized in that, The system includes the following modules: Module M1: Determine the starting point T of the time range of the batch observation data currently being processed s and the ending point T e ; Module M2: Recursively derive the orbits of cataloging targets to the corresponding time series The orbital positions and velocities in the J2000.0 inertial coordinate system, and the result is recorded as the orbit matrix O; Module M3: Based on the time series Generate the position interpolation function f with the orbital matrix O using the numpy / scipy package; Module M4: Input the current observation data to be associated. Module M5: Obtain Observation Time Series Measured values of observation data in the J2000.0 inertial coordinate system and the station position matrix P; Module M6: From the interpolation function f and the observed time series Generate the position matrix Q of the catalog target corresponding to the observed time series ; Module M7: Calculate the predicted value of the observable quantity of the position matrix Q relative to the station position matrix P. Module M8: Calculate the difference between the predicted value and the measured value of the observation data. Module M9: Calculate the associated object and confidence level based on the difference. Module M10: For the current batch of observation data to be processed, repeatedly call Module M4 to Module M9 until the end.

7. The correlation system for space debris observation data and cataloged orbit data according to claim 6, wherein In the said Module M2: The said time series is T s to T e The time series evenly distributed at equal time intervals within the range is a one-dimensional vector: Catalog target corresponding time series The orbit matrix O of is a three-dimensional matrix with dimensions M×6×N, where the 0th dimension corresponds to the time series The length M of the 0th dimension is the length of the time series The length of the 1st dimension 6 corresponds to position and velocity, and the length of the 2nd dimension N corresponds to the number of catalog targets; In the module M5: If the length of the observed time series is K, then the size of the station position matrix in the J2000.0 inertial coordinate system is K×3.

8. The correlation system for space debris observation data and cataloged orbit data according to claim 6, wherein The calculation of the predicted value of the observable quantity in Module M7 includes calculating the relative position matrix, distance matrix, and relative pointing matrix. The calculation of the predicted value of the observable quantity includes the following modules: Module M7.1: Calculate the relative position matrix R: R = Q - P (Formula 2) Where, the size of the target position matrix Q is K×3×N, and the size of the station position matrix P is K×3. When subtracting the two, the broadcasting function of the numpy array needs to be used. During the calculation, P is expanded to K×3×1 for operation, and the size of the resulting relative position matrix R is K×3×N. Module M7.2: Calculate the distance matrix ρ. The distance matrix ρ is the norm calculated along the first dimension of the relative position matrix R; the size of the distance matrix ρ is K×N. Module M7.3: Calculate the relative pointing matrix u: u = R / ρ (Formula 3) Where, the size of the relative position matrix R is K×3×N, and the size of the distance matrix ρ is K×N. When dividing the two, the broadcasting function of the numpy array needs to be used. During the calculation, ρ is expanded to K×1×N for operation, and the size of the resulting relative pointing matrix u is K×3×N.

9. The correlation system for space debris observation data and cataloged orbit data according to claim 6, wherein The calculation of the difference between the predicted value and the measured value of the ranging observable quantity in Module M8 includes calculating the distance difference matrix and the average distance difference within the observation duration. The calculation of the difference between the predicted value and the measured value of the angle measurement observable quantity includes calculating the angle measurement inner product matrix, the angle measurement difference matrix, and the average angle measurement difference within the observation duration. The calculation of the difference between the predicted value and the measured value of the ranging observable quantity includes the following modules: Module M8.1: Calculate the distance difference matrix ρ d : Among them, the distance matrix ρ is a two-dimensional matrix with a size of K×N. is a one-dimensional vector with a length of K. During the calculation, is expanded to K×1 for calculation, and the distance difference matrix ρ is calculated. d is a two-dimensional matrix with a size of K×N. Module M8.2: Calculate the average distance difference within the observation duration For the result of calculating the mean of the 0th dimension of the distance difference matrix ρ d is a one-dimensional vector with a length of N; The calculation of the difference between the predicted value and the measured value of the angle measurement observable quantity includes the following modules: Module M8.3: Calculate the angular inner product matrix u d : u d = sum(u·u m ) (Formula 5) Among them, the relative pointing matrix \(u\) is a three-dimensional matrix with a size of \(K\times3\times N\), which is the measurement pointing matrix \(u\). m , with a size of \(K\times3\). In the calculation, \(u\) m is extended to \(K\times3\times1\) for element-wise multiplication calculation. The result of \(u\cdot u\) m is still a three-dimensional matrix with a size of \(K\times3\times N\); the sum function represents summation along the first dimension of the three-dimensional matrix. The angular inner product matrix \(u\) d is a two-dimensional matrix with a size of \(K\times N\). Module M8.4: Calculate the angular difference matrix a d : a d = arccos(u d ) (Equation 6) where the arccos function represents the arccosine calculation, and the angular difference matrix a d has the same size as the angular inner product matrix u d which is K×N; Module M8.5: Calculate the average angular difference within the observation duration For the result of calculating the mean of the 0th dimension of the distance difference matrix a d is a one-dimensional vector with a length of N.

10. The correlation system for space debris observation data and cataloged orbit data according to claim 6, wherein The associated object calculation system in Module M9 includes calculating the difference evaluation factor, and splitting the difference evaluation factor to obtain the two smallest values, and judging the associated object and confidence level based on these two values. If the observation data only has angle measurement data, the evaluation factor is equal to the average angle measurement difference. If the observation data includes ranging and angle measurement values, the evaluation factor is calculated by adjusting the ranging difference and the angle measurement difference through a scaling factor. The calculation system of the association confidence level is: if the minimum value of the difference evaluation factor is less than the threshold and the second smallest value is greater than the threshold, the association confidence level is equal to 1; if both the minimum value and the second smallest value are less than the threshold, the association confidence level is calculated according to the formula. The associated object calculation system includes the following modules: Module M9.1: Calculate the differential evaluation factor is a one-dimensional vector with a length of N; Module M9.2: Split the difference evaluation factor by partition to obtain the two smallest values. Module M9.3: Obtain the minimum value j of the difference evaluation factor min and the second minimum value j 2nd ; Module M9.4: If the minimum value is greater than the association threshold, there is no matching object for this observation data. If the minimum value is less than the threshold, the observed data is associated with the cataloged orbit corresponding to the minimum value; If the observed data only contains angular measurement data, the evaluation factor is equal to the average angular measurement error If the observed data includes distance measurement and angular measurement values, the evaluation factor can be calculated by the following formula: Since radar equipment usually has high ranging accuracy but relatively poor angle measurement accuracy, the contribution of the angle measurement data to the evaluation factor is adjusted by the proportionality factor α, where α ≥ 0. When α = 0, the evaluation factor is completely provided by the ranging data; If the minimum value j min is less than the threshold and the second minimum value j 2nd is greater than the threshold, the associated confidence level c = 1. If the minimum value j min and the second minimum value j 2nd are both less than the threshold, the calculation system for the associated confidence level c is:

Images