Space target initial orbit determination method based on Gibbs algorithm

Abstract

The invention belongs to the technical field of aerospace orbit determination, and provides a space target initial orbit determination method based on a Gibbs algorithm in order to improve initial orbit determination precision and reduce errors, which comprises the following steps: acquiring a continuous position vector from an observation arc section; fitting the position vector of the observation arc section by using a least square method, and performing orbit plane fitting to obtain an orbit plane; position vector projection; carrying out mean value processing on projection vectors in the data set, and smoothing noise and abnormity in the data; according to the method, the situation that an existing Gibbs algorithm is only used for cognition of scarce measurement data is broken through, a geometric orbit determination method Gibbs method originally based on a pure vector of three pieces of distance measurement information is combined with an estimation algorithm based on least square fitting, and the space target initial orbit is determined. Initial orbit determination optimization is carried out, the orbit determination precision is improved, and the position error is minimized.

Description

Technical Field

[0001] The invention belongs to the technical field of aerospace orbit determination. Background Art

[0002] As human beings further explore space, the number of spacecraft and unknown space debris in orbit is increasing, and space safety issues are becoming more and more serious. Predicting the orbit of space objects in orbit and improving their prediction accuracy are of great significance for taking emergency collision avoidance measures for spacecraft in orbit and improving space safety. At the same time, it is also an important part of space situational awareness (SSA). At present, according to statistics from NASA ODPO (Orbital Debris Program Office), as of October 2023, most of the space debris larger than 10 cm is concentrated in low Earth orbit, and the number has exceeded 20,000. The number of space debris in other medium Earth orbits and high Earth orbits is also close to 5,000. On December 3, 2021, Starlink satellites approached the Chinese space station twice. For safety reasons, the Chinese space station took emergency collision avoidance measures.

[0003] In this context, it is crucial to strengthen the monitoring and tracking capabilities of space objects in orbit, and the determination and prediction of space target orbits is a major part of this. The determination of space target orbits generally includes three steps: (1) pre-processing of observation data; (2) initial orbit determination using short arc segment observation data; (3) orbit improvement using long arc segment. The initial orbit determination in step (2) above generally adopts a two-body model, and uses appropriate methods to make initial guesses of the orbit based on the available observation data. This is often done without prior information, and the orbit representation method is mainly based on the six orbital numbers, which are six parameters used to describe the size, shape, and position of the space target orbit.

[0004] Existing initial orbit determination methods can be roughly divided into two categories:

[0005] One is the deterministic method of directly calculating the six elements of the orbit in a scarce observation data environment. The observation data is usually at most 3. For example, the classic Gibbs three-vector orbit determination method in inertial space is extended to the relative motion space, and a Gibbs algorithm under relative motion conditions is given, such as Figure 1 As shown, the space target is in the space target orbit t 1 ,t 2 ,t 3 (t 1 <t 2 <t 3 ) The position vector r at three moments 1 、r 2 and r 3It can be found that the Gibbs algorithm uses pure vector analysis to calculate the velocity vector v at any time. 1 、v 2 and v 3 , and then calculate the six orbital elements.

[0006] The three-vector coplanar relationship: r 2 =c 1 r 1 +c 3 r 3 ,

[0007] Kepler orbital dynamics equation: r·e=h 2 / μ-r, e is the eccentricity, h is the specific angular momentum, μ is the gravitational constant, μ = 398600km 3 / s 2 ,

[0008] Combining the above formulas, we can get:

[0009] Where D = (r 1 × 2 )+(r 2 × 3 )+(r 3 × 1 ), N=(r 2 × 3 )r 1 +(r 3 × 1 )r 2 +(r 1 × 2 )r 3 ,

[0010] N is the ascending node vector. According to the relationship between the position vector and the velocity vector:

[0011] We can get:

[0012] Where: S = r 1 (r 2 -r 3 )+r 2 (r 3 -r 1 )+r 3 (r 1 -r 2 ).

[0013] The velocity vector v at a certain moment can be calculated by the above formula. If the position vector information at that moment is known, the orbit information of the space target can be obtained according to the six-element orbit calculation formula. However, the Gibbs algorithm is based on the premise that the above three-vector coplanar relationship is established. This coplanar characteristic has strict requirements and will lead to large errors in the initial orbit determination in actual observations. At the same time, the number of observation data available for the Gibbs algorithm is limited, which will limit the accuracy of the initial orbit determination.

[0014] The second is an estimation method in a tightly observed data environment, which usually uses complete data measured within a certain time range to optimize the state vector, rather than directly calculating it through mathematical algebra or vector methods. Most of them are fitting algorithms, and their optimization results depend on the length of the observed arc, the density of the data volume, and the algorithm itself. For example, the least squares method finds the optimal model parameters by minimizing the sum of squares of the errors between the actual observations and the model predictions. Although it uses relatively tight observations within the observable range, the computational cost is relatively high. Summary of the invention

[0015] In order to solve the technical problems in the prior art that the Gibbs algorithm is only applicable to scarce measurement data, resulting in low initial orbit determination accuracy, and the coplanar characteristics lead to large errors, and the calculation cost of the estimation algorithm is very high, the present invention proposes a "space target initial orbit determination method based on Gibbs algorithm" which combines the two.

[0016] like Figure 2 As shown, the method for determining the initial orbit of a space target based on the Gibbs algorithm includes the following steps:

[0017] Step 1: Collect observation arc point cloud data: Obtain continuous position vectors from the observation arc. Each position vector {r 1 ,r 2 ,…,r n} represents the spatial position of the target;

[0018] Step 2: Use the least squares method to fit the position vector of the observed arc segment and perform orbital plane fitting to obtain the orbital plane equation, x i ,y i , z i is a series of n observation points (n≥3), i=1,..,n,

[0019] The orbital plane equation is z = ax + by + c, c = 0,

[0020] Therefore, the least squares error function is constructed as:

[0021]

[0022] Consider a and b as independent variables, and write them into matrix form after differentiation:

[0023]

[0024] Substitute the position vector coordinates of the n observation points into the above formula, fit the values ​​of a and b, and substitute them into the orbital plane equation;

[0025] Step 3, position vector projection: Project all position vectors in the observation arc segment onto the orbital plane fitted in step 2, so that the three position vectors at different times are coplanar, and obtain the vector {R 1 ,R 2 ,…,R n};

[0026] like Figure 3 As shown, r is the initial position vector, n is the normal vector of the fitting orbital plane, R is the vector of r projected onto the fitting orbital plane, and r' is the vector of r projected onto n;

[0027] R=r-r',

[0028]

[0029] n=(a,b,-1), get the projection vector R of all position vectors in the observation arc segment;

[0030] Step 4: From the projection vector {R 1 ,R 2 ,…,R n}Select three sets of data sets S 1 , S 2 and S 3 Each set of data is distributed in different time ranges of the observation arc and does not overlap. The projection vectors in the data set are averaged to smooth the noise and anomalies in the data.

[0031] The position vector dataset of the Kth observation arc is S k ,in k = 1, 2, 3, j represents the jth observation value in each observation arc segment;

[0032] i represents the middle observation time of each observation arc segment, m is half the length of each observation arc segment, 2m+1 is less than or equal to the number of observations in each observation arc segment, and k=1, 2, 3, i, m∈N, i≥m, we get R 1 , R 2 and R 3 ;

[0033] Step 5: Use R 1 , R2 and R 3 As the input of Gibbs algorithm, the velocity vector v at the corresponding time is obtained 1 、v 2 and v 3 , using R 1 , R 2 and R 3 and v 1 、v 2 and v 3 Calculate the six orbital numbers to determine the initial orbit of the space target.

[0034] Technical effect:

[0035] The present invention breaks the existing perception that the Gibbs algorithm is only used for scarce measurement data. It combines the Gibbs method, a pure vector geometric orbit determination method originally based on three ranging information, with an estimation algorithm based on least squares fitting to perform initial orbit determination optimization, thereby finding a method for determining the initial orbit of space targets that is more suitable for current measurement data.

[0036] The present invention uses all the measurement data of the observed arc segment, the data is relatively compact, the data volume is large, and the orbit determination accuracy can be improved; at the same time, in view of the technical problem that the existing Gibbs algorithm has strict requirements on the coplanar characteristics of the measurement data, which leads to large errors, the least squares fitting and the mean algorithm are used to improve the skewness of the position vector and minimize the position error, so that the data is more in line with the plane characteristics. The mean algorithm can smooth the noise and anomalies in the data, increase the robustness of the Gibbs algorithm, improve the initial orbit determination accuracy while saving the calculation cost, and expand the application field of the existing Gibbs algorithm.

[0037] In order to better illustrate the technical effect of the present invention, the initial orbit determination accuracy and robust performance of the present method are verified by comparing the orbit determination success rate, the semi-major axis error and the algorithm running time. Data points at all observation times are selected in the observation arc for fitting, and three observation arcs are selected, namely, the 100s observation arc, the 300s observation arc and the 600s observation arc. The step lengths of the above observation arcs are all 1s. The results show that in terms of orbit determination success rate, in the 300s observation arc, the success rate of the present method is 91.9%, while the success rate of the existing Gibbs algorithm is 86.9%, an increase of 5% in success rate; in the 600s observation arc, the success rate of the present method is 98.7%, while the success rate of the existing Gibbs algorithm is 97%, an increase of 1.7%; in terms of the semi-major axis error, in the 100s observation arc, the existing Gibbs algorithm has more points gathered within an error range greater than 1000km, in the 300s observation arc, more points gathered within a range greater than 200km, and in the 600s observation arc, the semi-major axis error of the present method is 1.6%. The differences are basically concentrated in a range of less than 100 km, while the existing Gibbs algorithm has more points concentrated in a range of more than 100 km. This is because with the increase of the observation arc, any observation error and system error will accumulate and propagate over time, resulting in an increase in the uncertainty of the measurement result. The present invention can have a certain degree of smoothing effect on the uncertainty of the data and increase the algorithm's better orbit determination performance; in terms of the algorithm running time, the present invention does not increase the running time due to the addition of steps 2 to 4, and the calculation is concentrated below 8 ms, because the processing done is basically linear processing, and does not involve high-order term operations and matrix operations. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 This is the principle diagram of the three position vectors of the existing Gibbs algorithm.

[0039] Figure 2 The figure is a flowchart of the overall method of the present invention.

[0040] Figure 3 Schematic diagram illustrating the projection relationship of the position vector in step 3 of the present invention. DETAILED DESCRIPTION

[0041] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by ordinary technicians in the field without creative work by adopting the embodiments of the present invention are within the scope of protection of the present invention.

[0042] The specific derivation process of step 2 is:

[0043] x i,y i , z i is a series of n observation points (n≥3), i=1,..,n,

[0044] The orbital plane equation is: z=ax+by+c. Since the fitted orbital plane must pass through the center of the earth, c=0.

[0045] Therefore, the least squares error function is constructed as:

[0046] When the error function value s reaches the minimum value, the plane formed by a and b is the best fitting plane. Considering a and b as independent variables, the following equation can be obtained by taking their derivatives:

[0047]

[0048] make

[0049]

[0050] After finishing, we can get:

[0051]

[0052] Written in matrix form:

[0053]

[0054] Substitute the position vector coordinates of the n observation points into the above formula, fit the values ​​of a and b, and substitute them into the orbital plane equation.

[0055] In step 4,

[0056] i represents the middle observation time of each observation arc segment, m is half the length of each observation arc segment, 2m+1 is less than or equal to the number of observations in each observation arc segment, and k=1, 2, 3, i, m∈N, i≥m.

[0057] From the above formula, we can get R k The smoothing effect depends on S k The linear characteristics of S and the size of m. When the data interval of the observation arc segment is smaller, the size of m has a greater influence on S k The less sensitive the linear characteristics of S are, on the contrary, the size of m is k Therefore, when the data interval is smaller, the value of m can be appropriately increased to make R k The smoothing effect is the best. On the contrary, the value of m can be appropriately reduced.

[0058] The input data of the present invention is the position vector of the space target in the J2000 coordinate system. All methods that can obtain the target position vector by calculation or observation equipment can use the method described in the present invention to determine the orbit. At the same time, the present invention is not limited to the specific details in the above-mentioned implementation mode. Within the technical concept of the present invention, the technical scheme of the present invention can be subjected to a variety of simple modifications, and these simple modifications all belong to the protection scope of the present invention.

Claims

1. A method for determining the initial orbit of a space target based on the Gibbs algorithm, characterized in that: The steps include: Step 1: Collect observation arc point cloud data: Obtain continuous position vectors from the observation arc. Each position vector {r1, r2, …, r n } represents the spatial position of the target; Step 2: Use the least squares method to fit the position vector of the observed arc segment and perform orbital plane fitting to obtain the orbital plane equation, x i ,y i , z i is a series of n observation points (n≥3), i=1,..,n, The orbital plane equation is z = ax + by + c, c = 0, The least squares error function is constructed as: Consider a and b as independent variables, and write them into matrix form after differentiation: Substitute the position vector coordinates of the n observation points into the above formula, fit the values ​​of a and b, and substitute them into the orbital plane equation; Step 3, position vector projection: Project all position vectors in the observation arc segment into the orbital plane fitted in step 2, so that the three position vectors at different times are coplanar, and obtain the vector {R1, R2, …, R n }; r is the initial position vector, n is the normal vector of the fitting orbital plane, R is the vector projected from r to the fitting orbital plane, r ' is the vector projected from r to n; R＝r-r'， n=(a,b,-1), get the projection vector R of all position vectors in the observation arc segment; Step 4: From the projection vector {R1, R2, …, R n }, three data sets S1, S2 and S3 are selected, each of which is distributed in different time ranges of the observation arc and does not overlap. The projection vectors in the data set are averaged to smooth the noise and anomalies in the data. The position vector dataset of the Kth observation arc is S k ,in k = 1, 2, 3, j represents the jth observation value in each observation arc segment; i represents the middle observation time of each observation arc segment, m is half the length of each observation arc segment, 2m+1 is less than or equal to the number of observations in each observation arc segment, and k=1, 2, 3, i, m∈N, i≥m, and we get R1, R2 and R3; Step 5: Use R1, R2 and R3 as the input of the Gibbs algorithm to obtain the velocity vectors v1, v2 and v3 at the corresponding time. Use R1, R2 and R3 as well as v1, v2 and v3 to calculate the six orbital elements to determine the initial orbit of the space target.

2. The method for determining the initial orbit of a space target based on the Gibbs algorithm according to claim 1, characterized in that: The specific derivation process of step 2 is: x i ,y i , z i is a series of n observation points (n≥3), i=1,..,n, The orbital plane equation is: z = ax + by + c. Since the fitted orbital plane must pass through the center of the earth, c = 0. Therefore, the least squares error function is constructed as: When the value s of the error function reaches the minimum value, the plane formed by a and b is the best fitting plane; taking a and b as independent variables and taking their derivatives, we can get the following equation: make After finishing, we can get: Written in matrix form: Substitute the position vector coordinates of the n observation points into the above formula, fit the values ​​of a and b, and substitute them into the orbital plane equation.

Images