Multi-satellite cooperative localization algorithm and proximity configuration optimization

Abstract

The invention relates to the field of non-cooperative target positioning, and discloses a multi-sight cooperative positioning and multi-tracking star approaching configuration optimization algorithm in a non-cooperative target rapid approaching process. Firstly, an orbital kinetic equation of a target is used as a state equation; then, multiple tracking satellites use a spaceborne camera to measure the sight angle of the target relative to the tracking satellites, and a measurement model is established; and finally, estimating orbit parameters of the target based on an extended Kalman filter (EKF) algorithm. On the basis, the berthing orbit height of each tracking satellite and the true anomaly during orbital transfer serve as parameters to be optimized, the root-mean-square error (RMSE) estimated by the target position and the intersection duration t serve as target functions, an improved NSGA-II algorithm is used for optimizing, and the berthing orbit height of each tracking satellite and the true anomaly during orbital transfer are determined; and the target positioning precision is improved. The double-sight-line cooperative positioning system provided by the invention is high in observability, and meets the requirement of accurate positioning in the approaching process of non-cooperative targets in different planes; the approximate configuration is optimized by combining the designed improved NSGA-II algorithm, and the two tracking satellites have the advantages of high positioning precision, short task time consumption and the like under the optimal configuration.

Description

Technical Field

[0001] The invention discloses a cooperative positioning algorithm and a configuration optimization algorithm in a process of fast approach of non-cooperative targets on different surfaces, and relates to the field of non-cooperative target positioning. Background Art

[0002] Microsatellites for tasks such as active debris removal and on-orbit servicing have gradually become an international research hotspot. In these tasks, accurate positioning of non-cooperative targets is a necessary prerequisite. Using optical sensors based on line-of-sight positioning can meet the mass and power limitations of microsatellites. However, the observability of single line-of-sight positioning is poor, and using multiple satellites to observe the same target can effectively improve the observability of the system. In addition, the configuration of multiple tracking satellites during the approach process directly affects the performance of target positioning estimation. However, most of the existing studies assume that the tracking satellite and the target are in the same orbital plane, establish a state model based on relative dynamics, deduce the system error propagation law based on the change of geometric configuration, and analyze the influencing factors.

[0003] However, in the tasks related to the rapid approach of non-cooperative targets in different planes, the tracking satellite and the target are not in the same orbital plane, the relative distance is far, and the relative speed is large. Therefore, the traditional relative kinematic equations are no longer applicable, and a new line-of-sight measurement and positioning scheme is needed. On the other hand, the geometric configuration of multiple tracking satellites and the target during the approach process directly affects the performance of target positioning estimation. Therefore, the present invention proposes a configuration optimization algorithm based on an improved SGA-II algorithm, which improves the target positioning accuracy during the rapid approach of different plane targets. Summary of the invention

[0004] Purpose of the invention: The purpose of the present invention is to propose a multi-line of sight collaborative positioning algorithm for the precise positioning requirements of non-cooperative targets in different planes during rapid approach. In addition, the configuration of multiple tracking stars during the approach process directly affects the performance of target positioning estimation. Therefore, the present invention also proposes a configuration optimization algorithm based on the improved NSGA-II algorithm to improve the target positioning accuracy.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] (1) The orbital dynamics equation of the target is taken as the state equation.

[0007] The orbital parameters of the target in the J2000 coordinate system are composed of position and velocity, expressed as:

[0008] X(t) = [r, v] T (1)

[0009] Where r = [x, y, z] T is the target location, The target speed.

[0010] Considering the non-spherical gravity of the earth, the corresponding state equation is:

[0011]

[0012] Where μ is the Earth's gravitational constant, R e is the radius of the earth, r is the distance from the center of the spacecraft to the earth, J2 is the J2 perturbation coefficient, and f x 、f y 、f z It is the sum of the effects of the higher-order terms of the Earth's non-spherical perturbations, solar and lunar perturbations, solar pressure perturbations, and atmospheric perturbations.

[0013] (2) Multiple tracking satellites use onboard cameras to measure the target’s line-of-sight angle relative to themselves and construct a measurement model.

[0014] The tracking satellite uses the onboard camera to measure the target's line of sight angle relative to the tracker itself, that is, the azimuth angle β i and pitch angle α i , assuming that the position state of the tracking star i is (x i ,y i , z i ), the position state of the target is (x t ,y t , z t ), the measurement equation can be expressed as:

[0015]

[0016] In the formula, and is the measurement noise, and n is the number of tracking stars.

[0017] (3) Based on the state equation in step (1) and the measurement equation in step (2), the extended Kalman filter (EKF) algorithm is used to estimate the orbit parameters of the target.

[0018] (4) Taking the parking orbit altitude and true anomaly angle of each tracking satellite during orbit change as the parameters to be optimized, and the root mean square error (RMSE) of the target position estimation and the rendezvous time t as the objective function, the improved NSGA-II algorithm is used to optimize and determine the parking orbit altitude and true anomaly angle of each tracking satellite during orbit change, so as to improve the target positioning accuracy.

[0019]

[0020] Where f1 is the root mean square error of target position estimation, f2 is the intersection time, m is the number of target positioning times, and r is t is the true value of the target position, is the estimated value of the target position, and t is the intersection time.

[0021] (5) The NSGA-II algorithm has shortcomings such as poor search capability and slow convergence speed. Therefore, based on the characteristics of positioning performance optimization during the approach process, the algorithm is improved from three aspects:

[0022] The present invention uses the NDX operator to replace the traditional operator. Compared with the traditional operator, the NDX operator introduces a normally distributed random variable into the crossover link and performs a discrete recombination operation on the parent gene. The specific calculation formula is:

[0023]

[0024] In the formula, is the jth gene of the first offspring, is the jth gene of the second offspring, is the jth gene of the first parent, is the jth gene of the second parent, N(0,1) is a normally distributed random variable, r(j)∈(0,1) and obeys uniform distribution.

[0025] The present invention proposes that the variation distribution index changes in the manner of a cosine function, and the variation distribution index is in the interval [0, η max ] changes within the range, and gradually decreases with the increase of the number of iterations. The calculation formula is:

[0026]

[0027] Where η k is the variation distribution index of the k-th generation population, η max is the maximum mutation distribution index, k is the current iteration population size, and p is the maximum number of iterations.

[0028] When there are a large number of duplicate individuals in the population, the present invention uses a simulated disaster method to perform human intervention, by eliminating some duplicate individuals and reinitializing some population individuals, so that the algorithm continues to search for the global optimal solution.

[0029] The advantages of the present invention are:

[0030] 1. Using a satellite-borne camera to obtain the line-of-sight information of the target relative to the tracking star is a passive detection method. In addition, the dual-line-of-sight collaborative positioning system proposed by the present invention is not only highly observable, but also meets the positioning requirements for non-cooperative targets in different planes.

[0031] 2. The geometric relationship between the tracking satellite and the target during the approach process directly affects the target positioning performance. In the present invention, the two tracking satellites can start from different parking orbits, thereby changing the geometric configuration during the approach process, which lays the foundation for subsequent configuration optimization.

[0032] 3. Combined with the designed improved NSGA-II algorithm to optimize the approach configuration, the tracking satellite has the advantages of high positioning accuracy and short mission time under the optimal configuration. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] Figure 1 Flowchart for multi-view collaborative positioning

[0034] Figure 2 Optimize the flow chart for configuration;

[0035] Figure 3 This is a schematic diagram of a dual-line-of-sight collaborative positioning solution;

[0036] Figure 4 It is a schematic diagram of intersection;

[0037] Figure 5 is the target positioning error map; DETAILED DESCRIPTION

[0038] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention will be described in more detail below in conjunction with embodiments and drawings.

[0039] The present invention takes two tracking stars approaching the target star quickly as an example. Figure 3 This is a schematic diagram of the multi-line-of-sight collaborative positioning solution. The tracking satellite and the target are located in different orbital planes.

[0040] The detailed process includes: after the rendezvous timing is determined, the two tracking satellites maneuver from different parking orbits to quickly approach the target. During the approach process, the tracking satellite uses the onboard GNSS to obtain its own high-precision position information, and uses the onboard camera to measure the target's line of sight angle relative to itself, that is, the azimuth angle β 1 , β 2 and pitch angle α 1 , α 2 Combining the state equation with the measurement equation, the position and velocity of the target in the inertial system are estimated based on the extended Kalman filter algorithm to achieve accurate positioning of the target.

[0041] 1. Dual-line-of-sight collaborative positioning algorithm

[0042] The multi-line of sight collaborative positioning algorithm proposed in this invention is as follows: Figure 1 As shown in the figure. First, a state model based on the orbital dynamics model is established, then the line of sight angle of the non-cooperative target relative to the tracking star is obtained as the observation quantity, and then a measurement model of the collaborative positioning system based on line of sight information is established. Finally, the orbital parameters of the non-cooperative target are calculated based on the extended Kalman filter algorithm.

[0043] 1.1 Equation of state

[0044] The orbital parameters of the target in the J2000 coordinate system are composed of position and velocity, expressed as:

[0045] X(t) = [r, v] T (1)

[0046] Where r = [x, y, z] T is the target location, The target speed.

[0047] Considering only the harmonic term J2 which has the most significant impact on the target state, the corresponding state equation is:

[0048]

[0049] Where μ is the Earth's gravitational constant, R e is the radius of the earth, r is the distance from the center of the spacecraft to the earth, J2 is the J2 perturbation coefficient, f x 、f y 、f z It is the sum of the effects of the higher-order terms of the Earth's non-spherical perturbations, solar and lunar perturbations, solar pressure perturbations, and atmospheric perturbations.

[0050] 1.2 Measurement model

[0051] The tracking satellite uses the onboard camera to measure the target's line of sight angle relative to the tracker itself, that is, the azimuth angle β i and pitch angle α i , assuming that the position state of the tracking star i is (x i ,y i , z i ), the position state of the target is (x t ,y t , z t ), the measurement equation can be expressed as:

[0052]

[0053] In the formula, and To measure noise.

[0054] 1.3 Filtering algorithm

[0055] Linearize equation (2) and use time T as the discrete step length to obtain

[0056] X k+1 = ФX k +ΓW k (4)

[0057] In the formula, Ф=I+AT, G=[0 3×3 ;I 3×3], A is the Jacobian matrix of the state equation.

[0058] Let the state vector r t =(x t ,y t , z t ), expand the equation (3) by Taylor, and ignore the higher-order terms. After linearization, the observation matrix can be obtained:

[0059]

[0060] The EKF algorithm is used to coordinate the positioning of non-cooperative targets:

[0061]

[0062] In the formula, and P k，k-1 are the state vector prediction value and its prediction error variance matrix, and P k They are the estimated value of the state vector and its estimated error variance matrix, Ф k，k-1 is the state transition matrix from time k-1 to time k, Q k-1 is the system noise covariance matrix at time k-1, Γ k-1 is the system noise matrix at time k-1, R k is the measurement noise covariance matrix at time k, K k is the Kalman gain matrix at time k, H k is the observation matrix at time k, Z k is the measured value of the target's sight angle relative to the tracker at time k, is the observation function and I is the identity matrix.

[0063] 2. Approaching configuration optimization

[0064] The parking orbit heights of the two tracking satellites and the true anomaly angle during orbit change are taken as the parameters to be optimized, and the minimum root mean square error (RMSE) of the target position estimation and the minimum rendezvous time t are taken as the objective functions. As shown in formula (7), the improved NSGA-II is used to optimize and determine the parking orbit heights of the two tracking satellites and the true anomaly angle during orbit change.

[0065]

[0066] Where f1 is the root mean square error of target position estimation, f2 is the intersection time, m is the number of target positioning times, and r is i is the true value of the target position, is the estimated value of the target position, and t is the intersection time.

[0067] 3. Collaborative positioning performance optimization based on improved NSGA-II algorithm

[0068] The NSGA-II algorithm has shortcomings such as poor search ability and slow convergence speed. Therefore, based on the characteristics of positioning performance optimization in the approach process, the algorithm is improved from three aspects. The improved NSGA-II algorithm improves the configuration optimization process such as Figure 2 .

[0069] 3.1 NDX Operator

[0070] Compared with traditional operators, on the one hand, the normal distribution random variables are introduced into the crossover link, which expands the movement space of the offspring and the search space is wider, and is no longer limited to searching near the parent generation; on the other hand, the operation of discrete recombination of the parent generation's genes gives the new individual a wider value space, making it easier to jump out of the local optimum, further expanding the search space of the operator.

[0071]

[0072] In the formula, is the jth gene of the first offspring, is the jth gene of the second offspring, is the jth gene of the first parent, is the jth gene of the second parent, N(0,1) is a normally distributed random variable, r(j)∈(0,1) and obeys uniform distribution.

[0073] 3.2 Variable mutation distribution index

[0074] The variation distribution index in the traditional NSGA-II algorithm is a constant and cannot be changed appropriately with the evolution of the population. The present invention proposes that the variation distribution index changes in the manner of a cosine function, as shown in formula (9). The variation distribution index is in the interval [0, η max ] changes, and gradually decreases with the increase of the number of iterations, which improves the disadvantage of slow convergence caused by the wide search space of the NDX operator. It does not affect the efficiency of the global search in the early stage of the algorithm, but also speeds up the convergence in the later stage, making the Pareto solution set quickly approach the real frontier.

[0075]

[0076] Where η k is the mutation distribution index of the k-th generation population, η max is the maximum mutation distribution index, k is the current iteration population size, and p is the maximum number of iterations.

[0077] 3.3 Simulating disasters

[0078] When the NSGA-II algorithm generates offspring through crossover mutation, it generally generates a large number of duplicate individuals. Once the duplicate individuals have an advantage in the population, most offspring will retain the genetic characteristics of these duplicate individuals, which will cause the algorithm to mature prematurely and reduce the distribution diversity of the solution, which is not conducive to improving the optimization effect of the algorithm. In this regard, when there are a large number of duplicate individuals in the population, the present invention uses a simulated catastrophe method for human intervention. By eliminating some duplicate individuals and reinitializing some population individuals, the algorithm can jump out of the current solution space and continue to search for the global optimal solution.

[0079] 4. Simulation Verification

[0080] In order to verify the effectiveness of multi-line-of-sight collaborative positioning and multi-tracker approach configuration optimization algorithm in the process of fast approach of non-cooperative targets in different planes, digital simulation software was written in Matlab environment.

[0081] 4.1 Approaching configuration design

[0082] In order to reduce the energy consumption of the tracking satellite's orbit change, it is best not to perform maneuvers outside the orbital plane. Therefore, the intersection point is designed as the intersection of the non-cooperative target's orbit and the tracking satellite's orbital plane. In addition, the two tracking satellites have different initial parking orbits, but the orbital inclination and ascending node right ascension are the same. They change orbits at the same time to quickly approach the target. The schematic diagram is shown in the figure below. Figure 4 shown.

[0083] 4.2 Improved NSGA-II algorithm for configuration optimization

[0084] The simulation step is 1s. The optical sensor measurement accuracy: azimuth and pitch angle are both 0.0057° (1σ). The target is located in a circular orbit with an orbital inclination of 53°, a right ascension of the ascending node of 0°, and an orbital altitude of 1110km. The orbital inclination of the tracking satellite parking orbit is 58°, and the right ascension of the ascending node is 0°. The initial position error of the filter is 1000m (1σ) in all three axes; the velocity error is 10m / s (1σ) in all three axes.

[0085] The improved SGA-II algorithm is used to optimize the configuration, and the optimization results are shown in Table 1. The third group of configurations takes the shortest time, but the RMSE of the target position estimation is large; the RMSE of the configurations of other groups are similar, but the first group takes the shortest time. Therefore, the first group of configurations is selected as the optimal configuration.

[0086] Table 1 Configuration optimization results

[0087]

[0088] The simulation results of the target positioning three-axis error in both optimized and non-optimized configurations are shown in Figure 2. Figure 5As shown in the figure. The solid line in the figure is the target positioning error when the configuration is not optimized, and the dotted line is the target positioning error under the optimal configuration. Analysis shows that after the configuration is optimized, the dual-line of sight collaborative positioning algorithm converges after filtering for 600s, and the positioning error of the non-cooperative target is within 10m, and the position mean square error is 6.4m. Compared with the unoptimized one, the positioning accuracy is significantly improved and converges quickly, indicating that optimizing the approach configuration of the tracking satellite can improve the performance of positioning estimation.

[0089] In summary, the dual-line-of-sight collaborative positioning system proposed in the present invention has strong observability and meets the demand for precise positioning during the approach process of non-cooperative targets in different planes. Combined with the designed improved NSGA-II algorithm to optimize the approach configuration, experiments show that the two tracking satellites have the advantages of high positioning accuracy and short task time under the optimal configuration.

[0090] The above is only one embodiment of the present invention, which is only illustrative and not restrictive. Those skilled in the art understand that many changes, modifications, and even equivalents may be made to the present invention within the spirit and scope defined by the claims of the present invention, but all of them will fall within the scope of protection of the present invention.

Claims

1. A collaborative positioning and multi-tracker approach configuration optimization algorithm for non-cooperative targets in different planes during rapid approach, characterized by The following steps are involved: (1) The orbital dynamics equation of the target is used as the state equation; (2) Multiple tracking satellites use onboard cameras to measure the sight angle of the target relative to themselves and construct measurement equations; (3) combining the state equation in step (1) and the measurement equation in step (2), and using the extended Kalman filter (EKF) algorithm to estimate the orbit parameters of the target; (4) Taking the parking orbit height and the true anomaly angle of each tracking satellite during orbit change as the parameters to be optimized, and the root mean square error of the target position estimation and the rendezvous time as the objective function, the improved NSGA-II algorithm is used to optimize and determine the parking orbit height and the true anomaly angle of each tracking satellite during orbit change, so as to improve the target positioning accuracy.

2. A collaborative positioning method according to claim 1, characterized in that: The orbital parameters of the target in the J2000 coordinate system in step (1) are composed of position and velocity, expressed as: X(t)=[r,v] T (1) Where r = [x, y, z] T is the target location, The target speed. Considering the non-spherical gravity of the earth, the corresponding state equation is: Where μ is the Earth's gravitational constant, R e is the radius of the earth, r is the distance from the center of the spacecraft to the earth, J2 is the J2 perturbation coefficient, f x 、f y 、f z It is the sum of the effects of the higher-order terms of the Earth's non-spherical perturbations, solar and lunar perturbations, solar pressure perturbations, and atmospheric perturbations.

3. A collaborative positioning method according to claim 1, characterized in that: In step (2), the tracking star i uses the onboard camera to measure the line of sight angle of the target relative to the tracker itself, that is, the azimuth angle β i and pitch angle α i , assuming that the position state of the tracking star i is (x i ,y i , z i ), the position state of the target is (x t ,y t , z t ), the measurement equation can be expressed as: In the formula, and is the measurement noise, and n is the number of tracking stars.

4. A collaborative positioning method according to claim 1, characterized in that: Step (3) estimates the orbit parameters of the target using an extended Kalman filter (EKF) algorithm based on the state equation in step (1) and the measurement equation in step (2).

5. The multi-tracker approach configuration optimization algorithm according to claim 1, characterized in that: In step (4), the parking orbit height and true anomaly angle of each tracking satellite during orbit change are taken as the parameters to be optimized, and the root mean square error (RMSE) of the target position estimation and the rendezvous time t are taken as the objective function. The improved NSGA-II is used to optimize and determine the parking orbit height and true anomaly angle of each tracking satellite during orbit change. Where f1 is the root mean square error of target position estimation, f2 is the intersection time, m is the number of target positioning times, and r is i is the true value of the target position, is the estimated value of the target position, and t is the intersection time.

6. The improved NSGA-II algorithm according to claim 5, characterized in that: The NDX operator is used to replace the traditional operator, introduce the normal distribution random variable into the crossover link, and perform discrete recombination operations on the parental genes. The NDX operator calculation formula is: In the formula, is the jth gene of the first offspring, is the jth gene of the second offspring, is the jth gene of the first parent, is the jth gene of the second parent, N(0,1) is a normally distributed random variable, r(j)∈(0,1) and obeys uniform distribution.

7. The improved NSGA-II algorithm according to claim 5, characterized in that: In step (4), the variation distribution index changes in the manner of a cosine function. The present invention proposes that the variation distribution index is in the interval [0, η max ] changes in the manner of a cosine function, and gradually decreases with the increase in the number of iterations. The specific form is Where η k is the variation distribution index of the k-th generation population, η max is the maximum mutation distribution index, k is the current iteration population size, and p is the maximum number of iterations.

8. The improved NSGA-II algorithm according to claim 5, characterized in that: Step (4) adopts the method of simulating catastrophe. When there are a large number of duplicate individuals in the population, some of the duplicate individuals are removed, and some population individuals are reinitialized to continue searching for the global optimal solution.

Images