Target orbit determination method and system based on distributed filtering algorithm

Abstract

The present invention belongs to the field of aerospace technology, and specifically relates to a method and system for determining a target orbit based on a distributed filtering algorithm. The method comprises the following steps: S1, establishing a communication topology of a sensor network; S2, using a nonlinear dynamic model to represent the motion trajectory of a spacecraft and the position relationship between the spacecraft and a non-cooperative target; S3, using a measuring device to measure angle information and infer the position information of the non-cooperative target; S4, using a distributed filtering algorithm to filter and fuse the obtained non-cooperative target position information to obtain non-cooperative target state information; S5, using a distributed information fusion algorithm to perform distributed information fusion on the obtained non-cooperative target state information to determine the motion orbit of the non-cooperative target. The present invention has the characteristics of simple scheme, easy processing, small amount of numerical calculation, high-precision determination of the target orbit in a short time, and extremely high stability.

Description

Technical Field

[0001] The present invention belongs to the field of aerospace technology, and in particular relates to a target orbit determination method and system based on a distributed filtering algorithm. Background Art

[0002] In recent decades, with the continuous launch of rockets, the universe is filled with a large number of satellites and space debris launched by various countries. At the same time, with the continuous updating and development of detection technology, satellites are more sensitive to extraterrestrial meteorites, debris and UFOs. These foreign satellites, debris and UFOs are collectively referred to as non-cooperative space targets. Since non-cooperative space targets cannot be directly controlled, their uncertainty directly interferes with the normal operation of spacecraft. Therefore, in response to the problem that non-cooperative space targets cannot be observed, scientists from various countries have proposed situational awareness technology. Space situational awareness is an activity to obtain space target information through detection, analyze the situation of space targets and their impact on humans. Its main task is to accurately detect and track space targets, determine the characteristics of space targets that may pose a threat to the aerospace system, and analyze the target characteristic data. If these complementary and non-similar measurement and control systems, measurement data and orbit determination models are integrated to form a multi-source information fusion system for space target orbit determination, it can solve and make up for the shortcomings of single sensor observation technology, learn from each other's strengths and weaknesses, make full use of the information of each subsystem, improve the orbit determination accuracy, expand the scope of use, and improve the reliability of the orbit determination system.

[0003] Common target state estimation methods include Kalman filtering, least squares algorithm, particle filtering and other recognition estimation algorithms. Kalman filtering algorithm is most commonly used in target recognition and is also a linear filtering algorithm. The spatial target state determination using filtering algorithm is called spatial target tracking. However, since most observation systems and space target dynamic models are nonlinear, most existing literature uses nonlinear filtering algorithms such as extended Kalman filtering. Using nonlinear filtering algorithms can often obtain more accurate estimation results. The unscented Kalman filter (UKF) method is proposed based on the model and does not require linear approximation. It is more suitable for systems with strong nonlinearity, but it is sensitive to the observability of the system.

[0004] However, a single filter cannot guarantee the accuracy of target recognition. Combining Kalman filtering with information fusion methods, multiple satellite sensors work independently, each recognizes the target, and the satellites communicate with each other to fuse and identify the data. Compared with a single system, the multi-sensor system processes data more accurately and can well guarantee target recognition and tracking. There are two commonly used fusion structures. One is centralized fusion, in which all measurement data are sent to the center for direct processing. The advantage is that the information loss is minimal. The prior art proposes a robust centralized Kalman filter algorithm with uncertain noise variance. The disadvantage is that the amount of data processed is too large, which will cause serious computational load. The other structure is distributed fusion, in which the sensor generates a global optimal or suboptimal state estimate by using specific information fusion rules. In the prior art, the data fusion algorithm base containing convex combinations is integrated into the local Kalman filter, and the sigma point Kalman filter algorithm is proposed. The algorithm uses only multi-angle measurement system sampling values, providing a way of distributed fusion. In the prior art, the ADMM algorithm is also used to fuse the distributed Kalman filter algorithm to solve the constrained optimization problem, but the algorithm is complex, the calculation is cumbersome, and there are too many constraints, and the stability of the algorithm cannot be guaranteed.

[0005] Therefore, it is very necessary to design a target orbit determination method and system based on a distributed filtering algorithm that is simple, easy to handle, has low numerical calculation amount, can determine the target orbit with high precision in a short time, and has extremely high stability. Summary of the invention

[0006] The present invention aims to overcome the problem in the prior art that non-cooperative targets in space cannot be accurately observed. The present invention provides a method and system for determining a target orbit based on a distributed filtering algorithm, which is simple, easy to handle, has a small amount of numerical calculations, can determine the target orbit with high precision in a short time, and has extremely high stability.

[0007] In order to achieve the above-mentioned object of the invention, the present invention adopts the following technical solutions:

[0008] The target track determination method based on the distributed filtering algorithm includes the following steps:

[0009] S1, establish the communication topology of the sensor network, which is used to represent the communication connectivity between each satellite and each satellite;

[0010] S2, a nonlinear dynamic model is used to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target;

[0011] S3, using a measuring device to measure angle information and deduce the position information of the non-cooperative target;

[0012] S4, using a distributed filtering algorithm to filter and fuse the obtained non-cooperative target position information, to obtain the non-cooperative target state information after combining the observation data information of each satellite;

[0013] S5, using a distributed information fusion algorithm, performs distributed information fusion on the non-cooperative target state information obtained in step S4 to determine the motion trajectory of the non-cooperative target.

[0014] Preferably, step S1 comprises the following steps:

[0015] S11, define the mathematical expression of the directed graph as G = (V, E);

[0016] Where V = {1, 2, ..., n} represents the set of observed satellites; n is the total number of observed satellites; E (k) ={(i,j)|i,j∈V} is the set of network communication links, which is represented by the sum of other satellites that transmit information to satellite i at time k. If satellite i is connected to satellite j, they are called neighbor nodes.

[0017] S12, Definition

[0018] in represents the set of remaining satellites transmitting data to satellite i, represents the sum of neighboring satellites that transmit data outward from satellite i, It refers to the connected network formed by satellite i and its adjacent satellites. If there is a communication link sequence that enables communication between two satellites, the corresponding two satellites are considered to be connected; if any two satellites in a multi-source satellite system are connected, the communication topology of the corresponding system is considered to be a connected graph;

[0019] S13, in the topology diagram, the weighted and normalized adjacency matrix used is defined as:

[0020]

[0021] When satellite k is a neighboring satellite of satellite i, the weight coefficient is d ik The weight coefficient of satellite i itself is expressed as express.

[0022] Preferably, step S2 includes the following steps:

[0023] S21, a nonlinear dynamic model is used to represent the state transition. The model used is as follows:

[0024] x k =f(x k-1 )+w k-1 (30)

[0025] Among them, x k represents the state quantity at time k, the f function represents the state transfer equation in the model, and w k-1 represents the state noise at time k-1;

[0026] The noise is Gaussian white noise and satisfies the relationship cov(w k ,w j )=E(w k w j T )=Qδ jk , Q is the non-negative definite variance matrix of the observation noise, δ jk represents the unit impulse response, at w k With w j When it is uncorrelated, its value can be approximated to 1;

[0027] S22, assume that the satellite motion is not affected by external forces, the non-cooperative target t moves around the earth in a near-circular orbit, and the satellites i, j and k fly around the non-cooperative target t. At the same time, each satellite and the non-cooperative target satisfy the cw equation without external interference;

[0028] S23, set the relative distance between satellites i, j and k and non-cooperative target t to a, O-xyz represents the coordinate system centered on the center of the earth; without considering the influence of external forces, the relative dynamic equation between the spacecraft and the non-cooperative target is satisfied. Taking the i-th satellite as an example, the relative dynamic model cw equation is expressed by the following formula:

[0029]

[0030] Where x i represents the x-axis position of the ith satellite, z i represents the z-axis position of the ith satellite, represents the x-axis velocity of the ith satellite, represents the y-axis velocity of the i-th satellite, represents the z-axis velocity of the ith satellite, represents the x-axis acceleration of the ith satellite, represents the y-axis acceleration of the i-th satellite, represents the z-axis acceleration of the ith satellite, ω is the angular velocity vector, a is the distance between the satellite orbit and the target, μ is the earth's gravitational constant; when the initial relative state satisfies the condition When , the relative motion is periodic, and the motion trajectory of the non-cooperative target relative to the satellite is nearly circular;

[0031] S24, select the relative position and relative velocity of the space non-cooperative target relative to point O as the state quantity, that is, x tRepresents the x-axis position of the target t, y t represents the y-axis position of the target t, z t represents the z-axis position of the target t, represents the x-axis velocity of the target t, represents the y-axis speed of the target t, Represents the z-axis velocity of the target t, and the orbital relative dynamics equation is expressed as follows:

[0032]

[0033] Formula (4) can be further expressed as:

[0034]

[0035] S25, discretize and linearize the nonlinear dynamic model. The model is discretized with sampling time T, the system noise is w, and the discretized model is expressed as:

[0036] X (k) =X (k-1) +f(X (k-1) )T+w (k-1) (34)

[0037] Where: X (k) , X (k-1) , w (k-1) are the state quantity and model error at time k and k-1 respectively.

[0038] Discretize and linearize the cw equation, and the formula is as follows:

[0039]

[0040]

[0041] Where I is the 3×3 identity matrix, ρ = [x, y, z],

[0042] Finally, the expression of F is obtained. F is the discretized state transfer matrix, which can be expressed as:

[0043]

[0044] Preferably, step S3 includes the following steps:

[0045] S31, the state of a space satellite relative to a non-cooperative target is defined as Each satellite observes the elevation and azimuth angles θ, calculate the non-cooperative target position and transmit the target information to other connected satellites, define d x , d y, d z are respectively expressed as the relative distance between the satellite and the non-cooperative target at each coordinate, d x =|x t -x i |,d y =|y t -y i |,d z =|z t -z i |; Elevation and azimuth θ can be expressed as:

[0046]

[0047]

[0048] S32, the observation equation of a single satellite sensor is expressed as:

[0049]

[0050] in, It is represented as the observation value of the i-th sensor at time k; It is expressed as the observation noise of the ith sensor at time k and satisfies the equation cov(v k ,v j )=E(v k v j T )=Rδ jk , R is the non-negative definite variance matrix of the observation noise;

[0051] S33, when a single satellite observes a value, the observation equation needs to be linearized and discretized, which is specifically expressed by the following equation:

[0052] Z k =H k x k +v k (41)

[0053] The observation matrix H k It can be expressed by (14):

[0054]

[0055] Preferably, step S4 includes the following steps:

[0056] S41, without considering external interference, assume that satellite i and non-cooperative target t satisfy the relative dynamic equation, where the dynamic equation and measurement equation are expressed as follows:

[0057] X it(k+1)=F(k)X it (k) (43)

[0058]

[0059] Among them, X it (k) is the relative state between satellite i and target t, F(k) represents the discretized state transfer matrix, w it (k) is the measurement noise generated between satellite i and non-cooperative target t;

[0060] S42, after satellite i measures the angle information, an extended Kalman filter algorithm is used to perform filtering analysis to obtain an estimated target position.

[0061] Preferably, step S5 comprises the following steps:

[0062] S51, since the target t flies around the main spacecraft m, the relative dynamic equation is satisfied between the two, so the real trajectory of the non-cooperative target is expressed by the relative dynamic equation; each satellite first observes through its own observer, and then performs a filtering algorithm on its own data to obtain the optimal estimate of the state of the non-cooperative target, and then performs distributed information fusion on multiple data;

[0063] S52, in the distributed information fusion algorithm, the kinematic equation uses the relative dynamic equation between each satellite and the non-cooperative target, and the measurement equation is the target data observed by the satellite. The specific formula is as follows:

[0064]

[0065] Where Z i (k) is the observation equation, represents the estimated value of the non-cooperative target obtained by satellite i;

[0066] S53, observing the non-cooperative target separately according to the process of step S52, and then fusing the obtained data.

[0067] Preferably, the extended Kalman filter algorithm in step S42 includes the following steps:

[0068] S421, the extended Kalman filter algorithm adopts a nonlinear Kalman filter algorithm, and the algorithm process is as follows:

[0069]

[0070]

[0071] K k+1 =P k+1,k H(X k+1,k ) T[H(X k+1,k ) k+1,k ·H(X k+1,k ) T +R k+1 ] -1 (48)

[0072]

[0073] P k+1 =(IK k+1 H(X k+1,k ))P k+1,k (50).

[0074] X k+1,k represents the Kalman prediction value at time k+1 predicted by time k, represents the estimated value, K k+1 represents the Kalman coefficient, P k+1,k is the covariance matrix at time k+1 predicted by time k, Q k is the system measurement noise.

[0075] Preferably, the distributed filtering algorithm in step S4 uses information weight and information matrix to calculate the topological average value; S is set to represent the information weight, λ represents the information matrix, and the information weight and information matrix in the centralized filtering algorithm are represented as follows:

[0076]

[0077]

[0078] The superscript “+” is defined as an update operation. By using the update operation, any variable x is updated to the value of x. + Indicated by: Introducing consensus observation variables at the i-th sensor, combined with numerical values, represented by parameters q and B, the information weight and information matrix in the distributed filtering algorithm are represented as follows:

[0079]

[0080]

[0081]

[0082]

[0083] d ijis the weight coefficient between satellite i and its connected satellite j, q and B represent the improved information weight and information matrix after combining adjacent information errors, u is the information value of the combined observation value for calculating the information weight, and I is the influence factor containing the observation matrix value for calculating the information matrix.

[0084] The present invention also provides a target track determination system based on a distributed filtering algorithm, comprising:

[0085] A communication topology building module is used to build the communication topology of the sensor network, indicating the communication connectivity relationship between each satellite and between each satellite;

[0086] A nonlinear dynamics model building module is used to adopt a nonlinear dynamics model to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target;

[0087] A position information calculation module is used to measure angle information using a measuring device and calculate the position information of a non-cooperative target;

[0088] A filtering module is used to filter and fuse the obtained non-cooperative target position information using a distributed filtering algorithm to obtain the non-cooperative target state information after combining the observation data information of each satellite;

[0089] The fusion module is used to adopt a distributed information fusion algorithm to perform distributed information fusion on the obtained non-cooperative target state information and determine the motion trajectory of the non-cooperative target.

[0090] Compared with the prior art, the present invention has the following beneficial effects: (1) The present invention fuses the data of the target measured by multi-source satellite sensors by adopting an improved distributed Kalman filtering algorithm and only measuring the angle to observe the target dynamics. Compared with the observation estimation value of a single satellite sensor, the information is fused by adopting the improved distributed filtering algorithm with a consensus strategy, and the estimated value obtained is more accurate and has a smaller error; (2) The present invention simulates the observation and perception process of multiple satellites on unknown targets in outer space, estimates the target state through a filtering algorithm, and compares it with the actual state of the target, confirming the feasibility and accuracy of the filtering algorithm based on multi-source information fusion, and has higher accuracy than single satellite observation. BRIEF DESCRIPTION OF THE DRAWINGS

[0091] Figure 1 A schematic diagram of the motion trajectory of a spacecraft and the position relationship between the spacecraft and a non-cooperative target provided by an embodiment of the present invention;

[0092] Figure 2 A schematic diagram of a target observed by a single satellite observer provided in an embodiment of the present invention;

[0093] Figure 3A communication topology diagram of three satellites and a non-cooperative target provided by an embodiment of the present invention;

[0094] Figure 4 A schematic diagram of relative position error data of single filtering and distributed filtering on the x-axis provided by an embodiment of the present invention;

[0095] Figure 5 A schematic diagram of relative position error data of single filtering and distributed filtering on the y-axis provided by an embodiment of the present invention;

[0096] Figure 6 A schematic diagram of relative position error data of single filtering and distributed filtering on the z-axis provided by an embodiment of the present invention;

[0097] Figure 7 A schematic diagram of relative speed error data on the x-axis for single filtering and distributed filtering provided by an embodiment of the present invention;

[0098] Figure 8 A schematic diagram of relative speed error data on the y-axis for single filtering and distributed filtering provided by an embodiment of the present invention;

[0099] Fig. 9 A schematic diagram of relative speed error data of single filtering and distributed filtering on the z-axis provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0100] In order to more clearly illustrate the embodiments of the present invention, the specific implementation methods of the present invention will be described below with reference to the accompanying drawings. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those of ordinary skill in the art, other accompanying drawings and other implementation methods can be obtained based on these accompanying drawings without creative work.

[0101] Example:

[0102] The present invention provides a target track determination method based on a distributed filtering algorithm, comprising the following steps:

[0103] S1, establish the communication topology of the sensor network, which is used to represent the communication connectivity between each satellite and each satellite;

[0104] S2, a nonlinear dynamic model is used to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target;

[0105] S3, using a measuring device to measure angle information and deduce the position information of the non-cooperative target;

[0106] S4, using a distributed filtering algorithm to filter and fuse the obtained non-cooperative target position information, to obtain the non-cooperative target state information after combining the observation data information of each satellite;

[0107] S5, using a distributed information fusion algorithm, performs distributed information fusion on the non-cooperative target state information obtained in step S4 to determine the motion trajectory of the non-cooperative target.

[0108] The present invention adopts a multi-sensor distributed filtering algorithm. In order to clearly represent each satellite and the communication connectivity relationship between them, the present invention uses graph theory to represent the communication topology of the sensor network. The mathematical expression of the directed graph is defined as G = (V, E), where V = {1, 2, ..., n} represents the set of observed satellites, n is the total number of observed satellites, and E is the number of satellites. (k) ={(i,j)|i,j∈V} is the set of network communication links, which is represented by the sum of the remaining satellites that can transmit information to satellite i at time k. If satellite i is connected to satellite j, they are called neighbor nodes. Definition in represents the set of remaining satellites that can transmit data to satellite i, represents the total number of neighboring satellites to which satellite i can transmit data. It refers to the connected network formed by satellite i and its adjacent satellites. If there is a sequence of communication links that allow two satellites to communicate, then the two satellites are called connected. If any two satellites in a multi-source satellite system are connected, then the communication topology of the system is called a connected graph.

[0109] In the topological structure diagram, the weighted and normalized adjacency matrix used in the present invention is defined as:

[0110]

[0111] In this formula, when satellite k is a neighboring satellite of satellite i, the weight coefficient can be expressed as d ik Indicates that the weight coefficient of satellite i itself can be used It can be seen that the weight coefficient is a positive number, and the weight coefficient will not be 0 only when there is connectivity between satellites.

[0112] The present invention adopts a nonlinear dynamic model to represent the state transition. The model adopted by the present invention is represented by formula (2):

[0113] x k =f(x k-1 )+w k-1 (58)

[0114] Among them, x k represents the state quantity at time k, the f function represents the state transfer equation in the model, and w k-1represents the state noise at time k-1, specifically Gaussian white noise, and satisfies the relationship cov(w k ,w j )=E(w k w j T )=Qδ jk , Q is the non-negative definite variance matrix δ of the observation noise jk represents the unit impulse response, at w k With w j When it is uncorrelated, its value can be approximated to 1.

[0115] Assuming that the satellite motion is not affected by external forces, the target t moves around the earth in a near-circular orbit, and the satellite's trajectory relative to the earth is circular, such as Figure 1 As shown in Figure 1, satellites i, j and k fly around target t, and each satellite and target satisfy the cw equation without external interference. The relative distance between each satellite i, j and k and the non-cooperative target t is a, O-xyz represents the coordinate system centered on the center of the earth, and the spacecraft motion trajectory and the position relationship between the spacecraft and the non-cooperative target can be expressed as follows: Figure 1 express.

[0116] Without considering the influence of external forces, the spacecraft and the target satisfy the relative dynamic equation. Taking the i-th satellite as an example, the cw equation of the relative dynamic model can be expressed by equation (3):

[0117]

[0118] Where x i represents the x-axis position of the ith satellite, z i represents the z-axis position of the ith satellite, represents the x-axis velocity of the ith satellite, represents the y-axis velocity of the i-th satellite, represents the z-axis velocity of the ith satellite, represents the x-axis acceleration of the ith satellite, represents the y-axis acceleration of the i-th satellite, represents the z-axis acceleration of the ith satellite, ω is the angular velocity vector, a is the distance between the satellite orbit and the center of the earth, μ is the earth's gravitational constant; when the initial relative state satisfies the condition When the satellite is moving, the relative motion is periodic, and the motion trajectory of the non-cooperative target relative to the satellite is nearly circular.

[0119] The relative position and relative velocity of the non-cooperative target in space relative to point O are selected as the state quantity, that is, x t Represents the x-axis position of the target t, y t represents the y-axis position of the target t, zt represents the z-axis position of the target t, represents the x-axis velocity of the target t, represents the y-axis speed of the target t, represents the z-axis velocity of the target t. The orbital relative dynamics equation can be expressed as follows:

[0120]

[0121] Formula (4) can be further expressed as:

[0122]

[0123] In practical applications, it is necessary to discretize and linearize the nonlinear dynamic model. The model is discretized with sampling time T and system noise w. The discretized model can be expressed as:

[0124] X (k) =X (k-1) +f(X (k-1) )T+w (k-1) (62)

[0125] Where: X (k) , X (k-1) , w (k-1) are the state quantity and model error at time k and k-1 respectively.

[0126] Discretize and linearize the cw equation, where I is the 3×3 unit matrix, ρ = [x, y, z],

[0127]

[0128]

[0129] Finally, the expression of F is obtained. F is the discretized state transfer matrix, which can be expressed as:

[0130]

[0131] In order to obtain the relative state information of non-cooperative targets, the present invention needs to use visible light cameras, radars and other measuring equipment installed on satellites to measure angle information, thereby calculating the target position information. The state of the space satellite relative to the target is defined as Each satellite uses an angle-only observation system to observe the pitch angle and azimuth angle θ, calculate the position of the non-cooperative target and transmit the target information to other connected satellites, d x , d y , d z Respectively expressed as the relative distance between the satellite and the target under each coordinate, dx =|x t -x i |,d y =|y t -y i |,d z =|z t -z i |, elevation angle and azimuth angle θ can be expressed as:

[0132]

[0133]

[0134] For a more intuitive representation, the schematic diagram of a single satellite observer observing a target can be as follows: Figure 2 express, It is expressed as azimuth, pitch angle is expressed as θ, x, y, z represent the relative distance on x-axis, y-axis and z-axis respectively, the target distance is calculated by angle information, and the target speed is observed by speed observer. The angle measurement relationship can be expressed as Figure 2 express.

[0135] The observation equation of a single satellite sensor can be expressed as:

[0136]

[0137] in, It is expressed as the observation value of the i-th sensor at time k, It is expressed as the observation noise of the ith sensor at time k, which satisfies the equation cov(v k ,v j )=E(v k v j T )=Rδ jk , R is the non-negative definite variance matrix of the observation noise, δ jk represents the unit impulse response, at v k With v j When there is no correlation, its value can be approximated to 1. Each satellite observes the angle information through the sensor, calculates the target position information, and processes the target information obtained by each satellite through the fusion algorithm to obtain a more accurate value.

[0138] When observing a single satellite value, the observation equation needs to be linearized and discretized, which can be expressed by the following equation:

[0139] Z k =H k x k +v k (69)

[0140] The observation matrix H k It can be expressed by (14):

[0141]

[0142] The individual observation equations and relative kinematic equations of each satellite are established. It is assumed that the satellite cluster and non-cooperative targets fly around the main spacecraft and satisfy the relative kinematic equations. The satellite cluster and the target are not affected by other external forces. The orbit relative to the main satellite is a near-circular orbit. Each satellite is at a different position in the same altitude orbit and observes the target state information from different angles. While ensuring that the satellites can communicate and transmit data, a distributed filtering algorithm is used to filter and fuse the obtained target data to obtain the target state information after combining the observation data information of each satellite.

[0143] Without considering external interference, satellite i and target t satisfy the relative dynamic equation, where the dynamic equation and measurement equation can be expressed as follows: it (k) is the relative state between satellite i and target t, Z it (k) is the observation value of satellite i on target t:

[0144] X it (k+1)=F(k)X it (k) (71)

[0145]

[0146] Among them, X it (k) is the relative state between satellite i and target t, F(k) represents the discretized state transfer matrix, w it (k) is the measurement noise generated between satellite i and non-cooperative target t;

[0147] After satellite i measures the angle information, the extended Kalman filter algorithm is used to perform filtering analysis to obtain the estimated position of the target.

[0148] Since the target t flies around the main spacecraft m, the relative dynamic equation is satisfied between the two, so the true trajectory of the non-cooperative target can be expressed by the relative dynamic equation. Each satellite first observes through its own observer, and then performs a filtering algorithm on its own data to obtain the optimal estimate of the target state. Multiple data are then subjected to distributed information fusion to obtain more accurate data.

[0149] In the distributed information fusion algorithm, the kinematic equation uses the relative dynamic equation between the satellite and the target, and the measurement equation is the target data observed by the satellite:

[0150]

[0151] Where Z i (k) is the observation equation, Represents the target estimation value obtained by satellite i. According to the above scheme, non-cooperative targets are observed separately, and then their respective data are fused. This distributed filtering method not only takes into account the data processed by each satellite separately, but also in the communication topology diagram, each node can more accurately calculate the target status information, providing a basis for further processing of the target.

[0152] When a single satellite is used to observe and track a target, a traditional nonlinear filtering algorithm is used. This method is simple and fast and can roughly estimate the position state information. The filtering algorithm uses a nonlinear Kalman filtering algorithm, and its algorithm process is as follows:

[0153]

[0154]

[0155] K k+1 =P k+1,k H(X k+1,k ) T [H(X k+1,k ) k+1,k ·H(X k+1,k ) T +R k+1 ] -1 (76)

[0156]

[0157] P k+1 =(IK k+1 H(X k+1,k ))P k+1,k (78)

[0158] In this algorithm, the main factors affecting the target estimation state include the target covariance matrix and the Kalman gain.

[0159] The present invention adopts an improved distributed Kalman filter algorithm for the target state recognition scheme, and solves the algorithm performance problem of the distributed Kalman filter under the condition of low communication / computation complexity and without assuming local observability by combining the dynamic average consensus algorithm into the DKF.

[0160] The algorithm uses information weight and information matrix to calculate the topological average value. S represents the information weight and y represents the information matrix. Combining them with Kalman filtering can effectively solve the inconsistency and asynchrony problems in distributed algorithms. Equations (23) and (24) are the representation methods of information weight and information matrix in centralized filtering algorithms.

[0161]

[0162]

[0163] The superscript "+" is defined as an update operation, by using which any variable x is updated to its value by x. + For example, x(k+1)=Fx(k) can be rewritten as x + =Fx. A consensus observation variable is introduced at the i-th sensor. By combining the information matrix difference of the i-th sensor and the remaining sensors connected to it at the previous moment, the updated information weight and information matrix are more accurate. The combined values ​​are represented by q and B. The information weight and information matrix in the distributed filtering algorithm can be expressed by equations (27) and (28).

[0164]

[0165]

[0166]

[0167]

[0168] d ij is the weight coefficient between satellite i and its connected satellite j, q and B represent the improved information weight and information matrix after combining adjacent information errors, u is the information value of the combined observation value for calculating the information weight, and I is the influence factor containing the observation matrix value for calculating the information matrix.

[0169] Taking the i-th satellite as an example, formula (29) represents the information weight after combining the information weight of the connected satellites and the difference compensation at the previous moment. Formula (30) is the information matrix, which includes the information of the remaining satellites and the compensation at the previous moment.

[0170]

[0171]

[0172] By incorporating the information matrix and information weight into the Kalman filter algorithm, the target state estimation value and covariance matrix can be obtained.

[0173] M i =(P i -1 +S i ) -1 (87)

[0174]

[0175]

[0176]

[0177] P i + =FM i F T (91)

[0178] ε i represents the estimated value of the target by the adjacent satellites of satellite i after information fusion, represents the target estimate obtained by satellite i after combining with adjacent satellites, is the target prediction value at the next moment, P i + Represents the covariance matrix at the next moment.

[0179] In order to express the algorithm more clearly, the algorithm flow chart can be represented by the following figure:

[0180] Definition 1: If a matrix satisfies is full rank, j represents the set of sensors connected to sensor i, including sensor i, then sensor i can be said to be locally observable.

[0181] Assumption 1: The matrices F, H, R, Q are time-invariant and (A, H) is observable.

[0182] Assumption 2: The observed objects are distributed detectable, i.e., all sensors can capture the target dynamics.

[0183] Theorem 1: The target is observed by n satellites, and the target state is an m-dimensional matrix. Assume that the non-cooperative target is globally observable. And all satellites process the data individually and then fuse the processed data through the distributed algorithm (31)-(35). The improved distributed algorithm can solve the target object observation problem. Then the algorithm does not diverge, is in a stable state, and must satisfy:

[0184]

[0185] Where E[] represents the mathematical expectation,

[0186] prove:

[0187] In the proof process, the difference between the estimated value and the true value is expressed as η i Indicates that i is the i-th satellite, and the difference of n satellite values ​​is represented by a matrix, that is,

[0188]

[0189] The target estimate of the i-th satellite can be expressed by equation (37), where

[0190] Combining equations (32) and (34) into equation (37), we can obtain the expression of the target state at the next moment.

[0191]

[0192] Formula (38) can be rewritten as formula (39), where represents the Kronecker product, I m Represents an m-dimensional unit matrix. Aggregate n distributed data into one matrix to facilitate subsequent calculations, where as well as

[0193]

[0194]

[0195] Furthermore, the matrix FA is stable as shown in the previous assumption. Therefore, we can conclude that there exists a submultiplication matrix norm ∥·∥ p , so that ||FA|| p ≤ρ(FA)+ω 1 <1,||A T F T || p ≤ρ(FA)+ω 2 <1, ρ() represents the spectral radius of the matrix, which refers to the maximum value of the matrix eigenvalue. Expressing equation (39) as mathematical expectation, we can get the following expression:

[0196]

[0197] because but

[0198] The relationship between the spectral radius and the norm of can be expressed by (41):

[0199]

[0200] The norm has the following properties, ||X|| 2 ≤||X|| 1 ||X|| ∞ ,||D|| 2 ≤||D|| 1 ||D|| ∞ =1, then (41) can be further expressed.

[0201]

[0202] The convergence of the algorithm can be proved by formula (42): Since the matrix spectral radius is always less than 1, the algorithm is convergent.

[0203] Next, we prove the expected steady-state property of the algorithm, where h is the number of time steps. After h iterations, equation (40) can be expressed as

[0204]

[0205] After rearranging the inductive formula, we can obtain:

[0206]

[0207]

[0208] E[∈(j)] can be further written as formula (46), where

[0209]

[0210] Substituting equation (46) into (45),

[0211]

[0212] Considering h>r * , The value is close to 0, and equation (44) can be written as:

[0213]

[0214] Furthermore, the above equation can be rewritten as:

[0215]

[0216] Depend on You can get:

[0217]

[0218] Further simplifying the formula, we can get:

[0219]

[0220]

[0221] Therefore, it can be proved that the expectation of the difference between the estimated value and the true value obtained by the algorithm is in a convergence state, and the algorithm is feasible, ensuring that the expectation of the estimated value is less than a specific constant and will not diverge.

[0222] The present invention also provides a target track determination system based on a distributed filtering algorithm, comprising:

[0223] A communication topology building module is used to build the communication topology of the sensor network, indicating the communication connectivity relationship between each satellite and between each satellite;

[0224] A nonlinear dynamics model building module is used to adopt a nonlinear dynamics model to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target;

[0225] A position information calculation module is used to measure angle information using a measuring device and calculate the position information of a non-cooperative target;

[0226] A filtering module is used to filter and fuse the obtained non-cooperative target position information using a distributed filtering algorithm to obtain the non-cooperative target state information after combining the observation data information of each satellite;

[0227] The fusion module is used to adopt a distributed information fusion algorithm to perform distributed information fusion on the obtained non-cooperative target state information and determine the motion trajectory of the non-cooperative target.

[0228] In addition, the present invention verifies the feasibility of the distributed algorithm for non-cooperative space target state estimation in the aviation field through simulation experiments. In the simulation, it is assumed that three satellites in different positions in the same orbit track the non-cooperative target and collect data. At the same time, the observation satellite and the target fly around the same satellite, and their relative motion trajectory can be represented by the simulation diagram. At the same time, a fully connected communication topology structure is used to share target data and perform fusion calculations. The communication topology diagram can be obtained by Figure 3 express.

[0229] Assume that all three observation satellites use angle measurement to observe the target azimuth and pitch angle information, and use the Kalman filter algorithm to estimate the target state. The observation time is set to 300s, and the observation time interval is set to 1s. After that, the estimated target position and speed information are transmitted to adjacent satellites according to the communication topology for fusion calculation.

[0230] Assume the initial relative state of the target is i = 1, 2, 3, where the initial state is relative to the main satellite and i represents the observed satellite. The initial transition covariance matrix is ​​P i 0|0 =diag([1,1,1,1,1,1]), the system noise variance matrix is ​​Q=diag(10 -5 ,10 -5 ,10 -5 ,10 -5 ,10-5 ,10 -5 ), the measurement noise variance matrix is ​​R i =diag(1 2 ,1 2 ,1 2 ,0.1 2 ,0.1 2 ,0.1 2 ), the fusion weight coefficient d between different observation satellites ij Differently, for the main satellite that receives data from other satellites, the fusion weight coefficient accounts for a larger proportion, while the coefficient of the connected satellite is relatively small.

[0231] The simulation results are as follows Figures 4 to 9 As shown in the figure, the simulation results show that under the same initial covariance error and the same initial relative state of the target, the error between the estimated value and the true value calculated by the distributed fusion algorithm is in a dynamic equilibrium state, and the convergence performance has obvious advantages compared with the error value of the single sensor filtering algorithm. EKF represents single filtering, and DKF represents multi-source information distributed filtering. In the relative error diagram of the x-axis, the single filtering error value fluctuates around 10 meters in the first 50 seconds, and then the error decreases, and the accuracy reaches about 3 meters. From 100s to 150s, the accuracy is further improved, but after 150s, the error becomes larger again. It can be seen that the orbit determination accuracy of a single satellite observation target is unstable and the error fluctuates greatly. On the contrary, the distributed filtering algorithm integrating multiple satellites has very good stability. Similarly, when analyzing the x-axis velocity error, although the single filtering method is more accurate than the multi-source filtering method, there is still the problem of error fluctuation, and the stability of the observed target state cannot be guaranteed.

[0232] In the comparative analysis of the relative position and speed error of the y-axis, it can be obtained by Figure 5 and Figure 8 It can be seen that although the error accuracy on the y-axis is better than that on the x-axis, the single satellite filter has a larger fluctuation range in terms of position error. In the first 50 seconds, the maximum error reached 2.5 meters. In the following time, although the error continued to decrease, it could not be guaranteed to float at a certain stable value. In the multi-satellite distributed algorithm, although the error was sometimes high, the error coefficient could be stabilized at about 0.5 meters. Of course, in terms of the relative speed on the y-axis, the speed error of the distributed algorithm was stable at about 0.1m / s, but it was not as high as the error accuracy of the single satellite filter.

[0233] When analyzing the z-axis relative target state, the single filter error value is larger, and the intuitive comparison chart is more obvious. In the first 50 seconds, the single satellite filter error is as high as 8 meters. In the 50-200 seconds, the error fluctuates in the range of 0-4 meters. The error accuracy of the distributed filter can be stabilized at 0.5 meters. Similarly, on the z-axis, the speed error of the distributed algorithm is not as high as the error accuracy of the single satellite filter.

[0234] The simulation experiment shows that the distributed Kalman filter algorithm using fusion data has a more stable error in the target position and velocity estimation compared to the extended Kalman filter algorithm using a single satellite. As time goes by, there will be no large fluctuations in the value. Especially for position information, the distributed algorithm is obviously better for the position estimation of non-cooperative targets. When the accuracy of a single algorithm is 4 meters, the accuracy of the distributed algorithm can reach 0.5 meters. Although, in terms of the accuracy of velocity and position analysis, the accuracy of the distributed algorithm of multi-satellite numerical fusion is 0.05m / s less than that of the single filter algorithm, the difference in velocity error between the two is small and has no effect on target estimation. In practical applications, the accuracy of non-cooperative target position estimation is more important. If the target position information can be accurately determined in a short time, then the non-cooperative target can be further processed. Distributed filtering orbit determination can better meet this condition.

[0235] In view of the problem that non-cooperative targets in space cannot be accurately observed, the present invention proposes an improved distributed test orbit determination algorithm based on relative navigation. Compared with the calculation under the geocentric coordinate system, the scheme of the present invention is simple and easy to handle, with a small amount of numerical calculation, and can determine the target orbit with high precision in a short time. The orbit determination accuracy of the algorithm of the present invention in terms of position trajectory is significantly better than that of single satellite filtering orbit determination. At the same time, the improved filtering algorithm has extremely high stability.

[0236] The above description is only a detailed description of the preferred embodiments and principles of the present invention. For ordinary technicians in this field, according to the ideas provided by the present invention, there will be changes in the specific implementation methods, and these changes should also be regarded as the protection scope of the present invention.

Claims

1. A target track determination method based on a distributed filtering algorithm, characterized in that: The steps include: S1, establish the communication topology of the sensor network, which is used to represent the communication connectivity between each satellite and each satellite; S2, a nonlinear dynamic model is used to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target; S3, using a measuring device to measure angle information and deduce the position information of the non-cooperative target; S4, using a distributed filtering algorithm to filter and fuse the obtained non-cooperative target position information, to obtain the non-cooperative target state information after combining the observation data information of each satellite; S5, using a distributed information fusion algorithm to perform distributed information fusion on the non-cooperative target state information obtained in step S4 to determine the motion trajectory of the non-cooperative target; Step S4 includes the following steps: S41, without considering external interference, assume that satellite i and non-cooperative target t satisfy the relative dynamic equation, where the dynamic equation and measurement equation are expressed as follows: X it (k+1)=F(k)X it (k) (1) Among them, X it (k) is the relative state between satellite i and non-cooperative target t at time k, F(k) represents the discretized state transfer matrix at time k, and w it (k) is the measurement noise generated between satellite i and non-cooperative target t at time k; S42, after the satellite i measures the angle information, an extended Kalman filter algorithm is used to perform filtering analysis to obtain an estimated target position; The distributed filtering algorithm described in step S4 uses information weight and information matrix to calculate the topological average value; S is set to represent the information weight, λ represents the information matrix, and the information weight and information matrix in the centralized filtering algorithm are represented as follows: The superscript "+" is defined as an update operation; the consensus observation variable is introduced at the i-th satellite, and the numerical value is expressed by parameters q and B. The information weight and information matrix in the distributed filtering algorithm are expressed as follows: q and B represent the information weight and information matrix after the improved information error combined with the adjacent information, u is the information value of the combined observation value for calculating the information weight, I i is the influence factor of satellite i used to calculate the information matrix containing the observation matrix value.

2. The target track determination method based on the distributed filtering algorithm according to claim 1 is characterized in that: Step S1 includes the following steps: S11, the mathematical expression of the directed graph is defined as G = (v, e); where v = 1, 2, ..., n represents the set of observed satellites; n is the total number of observed satellites; e = (i, j) | i, j∈v represents the set of network communication links, e (k) It is represented as the sum of other satellites that transmit information to satellite i at time k. If satellite i is connected to satellite j, they are called neighbor nodes. S12, Definition in represents the set of remaining satellites transmitting data to satellite i, represents the sum of neighboring satellites that transmit data outward from satellite i, It refers to the connected network formed by satellite i and its adjacent satellites. If there is a communication link sequence that enables communication between two satellites, the corresponding two satellites are considered to be connected. If any two satellites in a multi-source satellite system are connected, the communication topology of the corresponding system is considered to be a connected graph. S13, in the topology diagram, the weighted and normalized adjacency matrix used is defined as: When satellite j is a neighboring satellite of satellite i, the weight coefficient is d ij The weight coefficient of satellite i itself is expressed as express.

3. The target track determination method based on the distributed filtering algorithm according to claim 2 is characterized in that: Step S2 includes the following steps: S21, a nonlinear dynamic model is used to represent the state transition. The model used is as follows: X (k) =f(X (k-1) )+w (k-1) ; Among them, X (k) , X (k-1) are the state quantities at time k and time k-1 respectively; the f(·) function represents the state transfer equation in the model; w (k-1) represents the system noise at time k-1, which is Gaussian white noise and satisfies the relationship cov(w (k-1) ,w (k -1) )=E(w (k-1) w (k-1)T )=Qδ jk , Q is the non-negative definite variance matrix of the system noise; δ jk represents the unit impulse response; S22, assume that the satellite motion is not affected by external forces, the non-cooperative target t moves around the earth in a near-circular orbit, satellites i and j fly around the non-cooperative target t, and at the same time, each satellite and the non-cooperative target satisfy the cw equation without external interference; S23, set a to represent the distance between the satellite orbit and the target, O-xyz represents the coordinate system centered on the center of the earth; without considering the influence of external forces, the relative dynamic equation between the spacecraft and the non-cooperative target is satisfied. Taking satellite i as an example, the relative dynamic model cw equation is expressed by the following formula: Where x i represents the x-axis position of satellite i, z i represents the z-axis position of satellite i, represents the x-axis velocity of satellite i, represents the y-axis velocity of satellite i, represents the z-axis velocity of satellite i, represents the x-axis acceleration of satellite i, represents the y-axis acceleration of satellite i, represents the z-axis acceleration of satellite i; ω is the angular velocity vector, μ is the Earth's gravitational constant; when the initial relative state satisfies the condition When , the relative motion is periodic, and the motion trajectory of the non-cooperative target relative to the satellite is nearly circular; S24, select the relative position and relative speed of the non-cooperative target t relative to the virtual center point O as the state quantity, that is, x t represents the x-axis position of the non-cooperative target t, y t represents the y-axis position of the non-cooperative target t, z t represents the z-axis position of the non-cooperative target t, represents the x-axis velocity of the non-cooperative target t, represents the y-axis velocity of the non-cooperative target t, represents the z-axis velocity of the non-cooperative target t; the orbital relative dynamics equation is expressed as follows: Formula (11) can be further expressed as: S25, discretize and linearize the nonlinear dynamic model, discretize the model with the system sampling time T, and the discretized model is expressed as: X (k) =X (k-1) +f(X (k-1) )T+w (k-1) (13) Where T is the system sampling time; the cw equation is discretized and linearized, and the formula is as follows: Where I is the 3×3 identity matrix, ρ = [x, y, z], Finally, the expression of F is obtained, where F is the discretized state transfer matrix, expressed as:

4. The target track determination method based on the distributed filtering algorithm according to claim 3 is characterized in that: Step S3 includes the following steps: S31, the state of a space satellite relative to a non-cooperative target is defined as Each satellite observes the elevation and azimuth angles θ, calculate the non-cooperative target position and transmit the target information to other connected satellites, define d x ,d y ,d z are respectively expressed as the relative distance between the satellite and the non-cooperative target at each coordinate, d x =|x t -x i |,d y =|y t -y i |,d z =|z t -z i |; The functional relationship between the state quantity and the observation quantity can be obtained by the pitch angle and azimuth angle θ represents: S32, the observation equation of a single satellite sensor is expressed as: In the formula, It is expressed as the observation value of the i-th satellite at time k; is the observation noise of satellite i at time k, and satisfies the equation R is the non-negative definite variance matrix of the observation noise; S33, when a single satellite observes a value, the observation equation needs to be linearized and discretized, which is specifically expressed by the following equation: Z k =H k x k +v k (20) The observation matrix H k It is expressed by the following formula:

5. The target track determination method based on the distributed filtering algorithm according to claim 1 is characterized in that: Step S5 includes the following steps: S51, since the target t flies around the main spacecraft m, the relative dynamic equation is satisfied between the two, so the real trajectory of the non-cooperative target is expressed by the relative dynamic equation; each satellite first observes through its own observer, and then performs a filtering algorithm on its own data to obtain the optimal estimate of the state of the non-cooperative target, and then performs distributed information fusion on multiple data; S52, in the distributed information fusion algorithm, the kinematic equation uses the relative dynamic equation between each satellite and the non-cooperative target, and the measurement equation is the target data observed by the satellite. The specific formula is as follows: in Denotes the estimated value of the non-cooperative target obtained by satellite i, w i (k) represents the observation noise; S53, observing the non-cooperative target separately according to the process of step S52, and then fusing the obtained data.

6. The target track determination method based on the distributed filtering algorithm according to claim 1 is characterized in that: The extended Kalman filter algorithm in step S42 includes the following steps: S421, the extended Kalman filter algorithm adopts a nonlinear Kalman filter algorithm, and the algorithm process is as follows: K k+1 =P k+1,k H(X k+1,k ) T [H(X k+1,k )P k+1,k ·H(X k+1,k ) T +R k+1 ] -1 (25) P k+1 =(I-K k+1 H(X k+1,k ))P k+1,k (27); X k+1,k represents the Kalman prediction value at time k+1 predicted by time k, represents the estimated value, K k+1 represents the Kalman coefficient, P k+1,k is the covariance matrix at time k+1 predicted by time k, Q k is the system measurement noise.

7. A target track determination system based on a distributed filtering algorithm, used to implement the target track determination method based on a distributed filtering algorithm according to any one of claims 1 to 6, characterized in that: The target track determination system based on the distributed filtering algorithm includes: A communication topology building module is used to build the communication topology of the sensor network, indicating the communication connectivity relationship between each satellite and between each satellite; A nonlinear dynamics model building module is used to adopt a nonlinear dynamics model to represent the trajectory of the spacecraft and the position relationship between the spacecraft and the non-cooperative target; A position information calculation module is used to measure angle information using a measuring device and calculate the position information of a non-cooperative target; A filtering module is used to filter and fuse the obtained non-cooperative target position information using a distributed filtering algorithm to obtain the non-cooperative target state information after combining the observation data information of each satellite; The fusion module is used to adopt a distributed information fusion algorithm to perform distributed information fusion on the obtained non-cooperative target state information and determine the motion trajectory of the non-cooperative target.

Images