A method of cooperative positioning of UAV swarms with low-orbit satellites based on information geometry

Abstract

The present invention relates to the field of low-orbit satellite electronic navigation, and is a multi-source collaborative fusion positioning method for large-scale unmanned aerial vehicle clusters in complex environments, and in particular to a low-orbit satellite collaborative positioning method for unmanned aerial vehicle clusters based on information geometry. The method utilizes the characteristics of high power of low-orbit satellite navigation signals and fast satellite position changes, combined with geometric ranging information between unmanned aerial vehicles in a unmanned aerial vehicle cluster, as well as inertial navigation and ADS‑B navigation systems, to overcome the positioning mutation error caused by the urban canyon effect, effectively improve the positioning accuracy of unmanned aerial vehicle clusters, and has good real-time performance, which can meet the navigation and positioning needs of large-scale unmanned aerial vehicle clusters in complex urban environments.

Description

Technical Field

[0001] The present invention relates to the field of low-orbit satellite electronic navigation, and is a multi-source collaborative fusion positioning method for large-scale UAV clusters in complex environments, and in particular to a UAV cluster low-orbit satellite collaborative positioning method (UAV Cluster Fusion Cooperative Positioning Algorithm With LEO Satellite System) based on information geometry. Background Art

[0002] The orbital altitude of low-orbit satellites is only 500-1200km. Compared with the global satellite navigation system, the orbital altitude of low-orbit satellites is lower, the propagation loss of satellite signals is smaller, and the positioning accuracy and anti-interference ability are higher. Therefore, in the field of high-precision positioning, the positioning method based on low-orbit satellites is being favored by researchers. However, due to the high flight speed of low-orbit satellites and the small coverage area of ​​a single satellite, ground drones cannot maintain long-term communication with a single satellite, making the communication link composed of the low-orbit satellite system and the drone cluster in rapid dynamic changes, increasing the difficulty of obtaining pseudo-range information; in addition, the rapid movement of low-orbit satellites will also make it impossible for a single drone to effectively obtain more than 4 pseudo-range data in complex environments.

[0003] The current low-orbit satellite positioning method is mainly based on Doppler single-star positioning, that is, through the high moving speed and signal power of low-orbit satellites, the Doppler effect of electromagnetic waves is used to solve the position of the target. Since the accuracy of the results of Doppler positioning alone is limited, the commonly used method is to combine the inertial navigation system for auxiliary positioning and use data fusion algorithms such as Kalman filtering to fuse data between navigation sources. As an excellent linear filtering method, Kalman filtering is widely used in data fusion, but its limitations are also obvious. When facing the scenario of cluster positioning, Kalman filtering needs to process multi-source data from multiple nodes. A lot of calculations are required to calculate the Kalman gain and update the estimation and error covariance, resulting in poor real-time performance of the algorithm.

[0004] At the same time, in the process of cluster positioning, in order to make full use of the geometric information in the cluster network, combined with intra-cluster collaborative positioning, the positioning accuracy can be greatly improved. At present, the intra-cluster collaborative positioning algorithm is mainly based on the neural network algorithm. However, in the face of complex network structures, the computational complexity of the neural network algorithm is relatively high, and the real-time positioning is very poor, which cannot meet the needs of cluster positioning. In addition, the existing patent "CN116931031A UAV Cluster Information Geometry Fusion Positioning Method Based on Low-orbit Satellites" mainly proposes a fusion method based on information geometry theory. This method solves the problem of poor real-time fusion positioning and makes large-scale cluster fusion positioning possible. However, this method does not take into account the characteristics of small coverage area and fast movement speed of low-orbit satellites, and cannot cope with the situation where the number of visible low-orbit satellites of the drone is less than 4 and cannot be positioned.

[0005] Therefore, there is an urgent need for a UAV cluster collaborative positioning method with good real-time performance, high positioning accuracy and suitable for low-orbit satellite systems to achieve real-time and accurate positioning of UAV clusters in urban environments. Summary of the invention

[0006] In view of the above problems, the present invention proposes a low-orbit satellite collaborative positioning method for UAV clusters based on information geometry. By utilizing the characteristics of high power of low-orbit satellite navigation signals and fast satellite position changes, combined with the geometric ranging information between UAVs in the UAV cluster, as well as inertial navigation and ADS-B navigation systems, it can overcome the positioning mutation errors caused by the urban canyon effect, effectively improve the positioning accuracy of UAV clusters, and has good real-time performance, which can meet the navigation and positioning needs of large-scale UAV clusters in complex urban environments.

[0007] The technical solution of the present invention is:

[0008] The method for collaborative positioning of UAV cluster low-orbit satellites based on information geometry is characterized by comprising the following steps:

[0009] Step 1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory, ensure that the number of pseudo-ranges of drones in the cluster is 4, and solve the three-dimensional covariance matrix of the relative positioning error;

[0010] Step 2: According to the navigation principles and error transfer models of low-orbit satellites, inertial navigation and ADS-B navigation, solve the positioning results and their three-dimensional covariance matrices of low-orbit satellites, inertial navigation and ADS-B navigation;

[0011] Step 3: Compare the positioning results of the low-orbit satellite for the UAV A to be positioned obtained in step 1 and the collaborative UAV A obtained in step 2 i The relative positioning results of the target A are fused to obtain the fused positioning results, and a UAV cluster information geometric fusion positioning model is constructed.

[0012] Furthermore, the step 1 is specifically as follows:

[0013] Step 1-1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory: equip each drone in the cluster with a low-orbit satellite receiver, ranging equipment, inertial navigation, and ADS-B navigation to ensure that each navigation system can work independently. The drone can communicate with the low-orbit satellite and establish communication with other drones within the communication distance, and then measure the current distance and direction to the other party;

[0014] Step 1-2: Process of obtaining the positioning of low-orbit satellites: Each drone in the drone cluster communicates with the low-orbit satellite through its own low-orbit satellite receiver, and obtains the low-orbit satellite navigation message and pseudo-range data through pseudo-range multiplexing of surrounding drones, constructs a positioning solution equation group based on the low-orbit satellite pseudo-range, and uses the least squares method and Newton iteration method to solve the equation group to solve the low-orbit satellite positioning coordinates of each drone in the cluster;

[0015] Step 1-3: According to the positioning error introduced by relative positioning during pseudorange reuse, the three-dimensional covariance matrix of the low-orbit satellite positioning results of each UAV is obtained by analyzing the error transmission model of the nodes in the factor graph.

[0016] Optionally, any target node m to be located q As an example, an internal factor graph is established. Since the coordinate axes are independent of each other in the coordinate solution process of the node to be located, the three-dimensional positioning problem is transformed into a one-dimensional problem on three branches, which is divided into x-coordinate, y-coordinate and z-coordinate. The three parts are connected by Euler's formula and expressed as:

[0017]

[0018] The constraint relationship formula is:

[0019]

[0020] in,

[0021] Δx iq , Δy iq and Δz iq They represent the relative distances between the node to be located q and its neighbor node i on the coordinate axis respectively;

[0022] d iq Indicates the relative distance between nodes;

[0023] (x i ,y i ,z i ) represents the coordinate value of neighbor node i;

[0024] (xq ,y q ,z q ) represents node m q The coordinates of .

[0025] Furthermore, in step 2, the positioning result of the low-orbit satellite is:

[0026]

[0027] in,

[0028] x k , x k-1 are the calculation results of the kth and k-1th Newton iterations respectively;

[0029] Δx, Δy, Δz, Δτ are the differences between the two solutions in the three directions during the Newton iteration process and the differences between the two clock error solutions.

[0030] τ is the clock difference between the low-orbit satellite and the ground drone.

[0031] Furthermore, in step 2, it is assumed that after the inertial navigation correction, p n =(x,y,z) T At this moment, the positioning result of the drone to be positioned is p n =(x,y,z) T , then the information geometry probability model of the inertial navigation source is:

[0032] p INS ~N(p n -μ INS ,∑ INS )

[0033] in,

[0034]

[0035] in,

[0036] B a =(B ax ,B ay ,B az ) T is the static zero bias of the accelerometer;

[0037] B ω =(B ωx ,B ωy ,B ωz ) T is the static bias of the gyroscope;

[0038] is the standard deviation of the accelerometer noise;

[0039] is the standard deviation of the gyroscope noise;

[0040] 1 / T s is the sampling frequency;

[0041] [f n ×] is The antisymmetric matrix of the accelerometer measurements under the system, and Direction cosine matrix By G b Sure;

[0042] [f n ×].^2 indicates that [f n ×] Square each element of the matrix;

[0043] diag(·) is a diagonal matrix construction function.

[0044] Furthermore, in step 2, the information geometry probability model of the ADS-B navigation positioning source can be expressed as:

[0045] p RNS ~N(μ RNS ,∑ RNS )

[0046] in,

[0047] μ RNS is the location information of the node to be located;

[0048] ∑ RNS is the covariance matrix of the location information of the node to be located.

[0049] Furthermore, in step 3, the process of obtaining the fusion positioning result is:

[0050] The positioning results from low-orbit satellites, inertial navigation, and ADS-B navigation are constructed as three-dimensional Gaussian probability density functions. Combined with the Riemann manifold under the Gaussian distribution family, the KL divergence of the probability density function of the positioning results on the Riemann manifold is calculated to approximate its Fisher information distance. When the sum of the Fisher information distances from a point in the Riemann manifold to each information source is the shortest, the probability density function corresponding to the point is the fused positioning result, thereby achieving accurate positioning.

[0051] Optionally, in step 3, the fusion positioning result is:

[0052]

[0053] in,

[0054] ω i is the weight factor, 0≤ωi ≤1 and satisfy the constraints ω1+ω2+…+ω N =1;

[0055] μ i and∑ i Represents the mean and covariance matrix of the probability density distribution to be fused;

[0056] μ f ,∑ f Represents the mean and covariance matrix of the fusion result.

[0057] Beneficial Effects

[0058] 1. The present invention proposes a method for cooperative positioning of UAV clusters using low-orbit satellites based on information geometry. By taking advantage of the higher power of low-orbit satellite signals to the ground, the positioning mutation error caused by the urban canyon effect in complex environments can be significantly suppressed, and the problem of being unable to position due to insufficient number of low-orbit satellite pseudoranges can be solved, thereby improving the positioning accuracy of UAV clusters in urban environments.

[0059] 2. The present invention converts low-orbit satellite positioning information, inertial navigation and ADS-B navigation positioning information into a probability density model to achieve rapid fusion of multi-source positioning information, thereby significantly improving the real-time and stability of UAV cluster positioning;

[0060] 3. The present invention uses the communication capability between drone clusters and the ranging lateral information within the cluster to complete pseudo-range multiplexing between drones when the number of visible satellites is insufficient, so that drones can use their own pseudo-range data and the pseudo-range data of surrounding drones to complete pseudo-range positioning, solving the problem of being unable to position when the number of visible satellites is less than 4. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 Co-locating internal node factor graph model in the embodiment;

[0062] Figure 2 It is a model of low-orbit satellite collaborative positioning system of UAV cluster based on information geometry in the embodiment;

[0063] Figure 3 The real-time comparison between the method of the present invention and the traditional method in the embodiment;

[0064] Figure 4 The stability comparison between the method of the present invention and the traditional method in the face of mutation errors in the embodiment;

[0065] Figure 5 The figure is a comparison of the positioning accuracy of the method of the present invention and the traditional method for different cluster sizes in the embodiment. DETAILED DESCRIPTION

[0066] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, and should not be construed as limiting the present invention.

[0067] This embodiment adopts a UAV cluster low-orbit satellite collaborative positioning method based on information geometry proposed by the present invention to perform real-time and accurate positioning of a UAV cluster with n=10 UAVs in a scene size of 5000m×5000m×500m, and the communication distance between UAVs is L=2000m, including the following steps:

[0068] Step 1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory, ensure that the number of pseudo-ranges of drones in the cluster is 4, and solve the three-dimensional covariance matrix of the relative positioning error, as follows:

[0069] Step 1-1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory: equip each drone in the cluster with a low-orbit satellite receiver, ranging equipment, inertial navigation, and ADS-B navigation to ensure that each navigation system can work independently. The drone can communicate with the low-orbit satellite and establish communication with other drones within the communication distance, and then measure the current distance and direction to the other party;

[0070] Step 1-2: Process of obtaining the positioning of low-orbit satellites: Each drone in the drone cluster communicates with the low-orbit satellite through its own low-orbit satellite receiver, and obtains the low-orbit satellite navigation message and pseudo-range data through pseudo-range multiplexing of surrounding drones, constructs a positioning solution equation group based on the low-orbit satellite pseudo-range, and uses the least squares method and Newton iteration method to solve the equation group to solve the low-orbit satellite positioning coordinates of each drone in the cluster;

[0071] Step 1-3: According to the positioning error introduced by relative positioning during pseudorange reuse, the three-dimensional covariance matrix of the low-orbit satellite positioning results of each UAV is obtained by analyzing the error transmission model of the nodes in the factor graph.

[0072] like Figure 1 As shown, take any target node m to be located q As an example, an internal factor graph is established. In the collaborative positioning method based on the factor graph, first, the adjacent positioning targets are collected for the node m to be positioned. q coordinate estimation information of the node to be located m; secondly, each positioning target provides the function node with the coordinate estimation information of the node to be located m q The required reliability information; finally, the reliability information is iteratively updated on the factor graph until convergence.

[0073] Since the coordinate axes are independent of each other in the process of solving the coordinates of the node to be located, the three-dimensional positioning problem is transformed into a one-dimensional problem on three branches, which is divided into x-coordinate, y-coordinate and z-coordinate. The three parts are connected by Euler's formula and expressed as:

[0074]

[0075] Where Δx iq , Δy iq and Δz iq They represent the relative distances between the node to be located q and its neighbor node i on the coordinate axis; d iq Indicates the relative distance between nodes.

[0076] The constraint relationship formula is:

[0077]

[0078] Among them, (x i ,y i ,z i ) represents the coordinate value of neighbor node i, (x q ,y q ,z q ) represents node m q The coordinates of .

[0079] In the pseudo-range migration algorithm based on factor graph, "credibility information" is iteratively transmitted between variable nodes and function nodes. Here, credibility information is mainly expressed as mean and variance information. Taking the kth iteration as an example, the information transmission between variable nodes and function nodes is as follows:

[0080] The ranging information is provided by the function node D i Enter the internal factor graph, represented by node D i The generated variable d i , subject to the mean The variance is The Gaussian distribution of

[0081]

[0082] After the ranging information enters the factor graph, the function node E i Convert q x Coordinates, q y Coordinates and q z Confidence information between coordinates. Assume that the surface equation is F(Δx, Δy, Δz, d) = Δx 2 +Δy 2 +Δz 2 -d 2 , at point (Δx iq ,Δy iq,Δz iq ,d iq ), the tangent equation of the surface can be simplified to:

[0083] Δx iq Δx+Δy iq Δy+Δz iq Δz-d iq d=0……………………(4)

[0084] The above formula is the linear constraint relationship between Δx, Δy, Δz, and d. Taking Δx as an example, the variable node Δx obeys a Gaussian distribution with a mean of m. Δx , then the variance of Δx

[0085]

[0086] Confidence Information are all Gaussian, and the confidence information from the function node to the variable can be expressed as:

[0087]

[0088] In obtaining Δx iq After the credibility information is obtained, the function node A is used i Realize the conversion between relative position and absolute position information, as shown in the following formula:

[0089]

[0090]

[0091] in, They are from the variable nodes The variance of the Gaussian reliability information.

[0092] After the upward iteration is completed, the downward iteration is used to prepare for the next update, and the variable node x q Based on the calculation criteria of the sum-product algorithm, update the function node A i The reliability information of

[0093]

[0094] Combining the characteristics of multiple Gaussian distribution variables that still obey Gaussian distribution after multiplication, the positioning of the target node to be located in the collaborative network is realized, and the variable node To the function node Taking confidence update as an example, the update process is as follows:

[0095]

[0096]

[0097]

[0098] Step 2: According to the navigation principles and error transfer models of low-orbit satellites, inertial navigation and ADS-B navigation, solve the positioning results and their three-dimensional covariance matrices of low-orbit satellites, inertial navigation and ADS-B navigation;

[0099] like Figure 2 As shown in the figure, it is assumed that UAV A1 has a total of n visible low-orbit satellites, and the signal-to-noise ratio of each low-orbit satellite signal is SN i (i=1,…,n), the low-orbit satellite signal with the largest signal-to-noise ratio is defined as the optimal, and the optimal pseudorange selection criterion is:

[0100] [min_SN,i]=max(SN i )………………………………(13)

[0101] Among them, min_SN is the signal-to-noise ratio of the optimal low-orbit satellite, i is the serial number of the low-orbit satellite, and the optimal pseudorange of UAV A1 comes from the low-orbit satellite L i .

[0102] Assuming that the number of drones that drone A1 can communicate with is m, the optimal pseudorange of the communicative drone can be transferred laterally through the ranging between drones, and drone A1 can indirectly use the pseudorange data of the communicative drone. Therefore, drone A1 can use a total of m+1 pseudorange data, and all pseudorange data are constructed into the pseudorange data matrix of drone A1. Through the minimum GDOP criterion of pseudorange, the four matrices with the best geometric properties and the most suitable for drone A1 positioning are selected from the pseudorange data matrix.

[0103] Assume that the coordinates of the four low-orbit satellites are S i =(X i ,Y i ,Z i )(i=1,2,3,4), low orbit star S i With drone A i The pseudo-range information is ρ i (i=1,2,3,4), and drone A i The real position coordinates of (i=1,2,3) are P i =(x i ,y i ,z i ) T , the coordinates of UAV A4 are P = (x, y, z) T , the propagation speed of electromagnetic waves in vacuum is a constant c, and the clock difference between the low-orbit satellite and the ground drone is τ, that is,

[0104]

[0105] Substitute the relative positioning result in step 1 to obtain:

[0106]

[0107] After the pseudo-range shift of the factor graph, Δx iq , Δy iq , Δz iq The variance will be converted into the position ambiguity of X1, Y1, and Z1, which will in turn affect the positioning results of the drone cluster.

[0108] Rearrange the above formula and let

[0109]

[0110] Linearize the above equation and get

[0111] G[Δx Δy Δz c·△τ] T =b……………………(17)

[0112] in,

[0113]

[0114]

[0115] I (x)4 (x k-1 ) is r i The partial derivative in the x direction is k-1 Negative values ​​of x k-1 is the calculation result of the k-1th Newton iteration, Δx, Δy, Δz, Δτ are the differences in three directions between the two solutions in the Newton iteration process and the differences between the two clock error solutions. It can be seen that the Jacobian matrix G is only related to the positional relationship between the low-orbit satellite and the drone.

[0116] Using the least squares method, we can get

[0117] [Δx Δy Δz c·△τ] T =(G T G) -1 G T b…………………………(20)

[0118] Then the positioning result after this Newton iteration is

[0119]

[0120] The drone is equipped with an inertial navigation source. The output of the inertial sensor accelerometer and gyroscope can be approximately represented by the following mathematical model:

[0121]

[0122] Among them, the superscript b represents the drone coordinate system, and the subscripts a and ω are the relevant quantities of the accelerometer and gyroscope, respectively. and represents the output value, f and ω represent the true value, S represents the proportional factor error, N represents the non-coincidence and non-orthogonality error, B represents the zero bias, and ε represents the random noise.

[0123] Assume that the drone is stationary, the accelerometer and gyroscope noises are independent of each other, and the accelerometer bias and attitude error are independent of each other, that is, the accelerometer measurement value is:

[0124] f b =[0 0-g] T …………………………………(twenty three)

[0125] When a static bias is present in the accelerometer, it introduces an error in velocity that is proportional to time, and an error in position that is proportional to the square of time:

[0126]

[0127] When there is a static zero bias in the gyroscope, a yaw error proportional to time is introduced:

[0128]

[0129] At the same time, after solving the INS update equation, it will cause speed and position errors in the horizontal direction:

[0130]

[0131] Therefore, the gyro constant error will introduce a constant in the velocity that is proportional to t 2 and introduces an error in the horizontal position proportional to t 3 The velocity measurement error caused by the integration of the accelerometer white noise is

[0132]

[0133] in, is the standard deviation of the accelerometer noise.

[0134] The white noise of the accelerometer will cause a second-order random walk error in the position measurement, with a mean of 0 and a standard deviation of:

[0135]

[0136] After integrating the gyroscope white noise, we can get the zero-mean angle random walk error, and its standard deviation is:

[0137]

[0138] in, is the standard deviation of the gyroscope noise.

[0139] The white noise of the gyroscope will cause a second-order random walk error in the horizontal velocity measurement and a third-order random walk error in the horizontal position measurement. The mean of the error is 0 and the standard deviation is:

[0140]

[0141] Assume that at time t after the inertial navigation correction, the attitude of the UAV to be positioned is G b =(φ,θ,ψ) T , the output of the accelerometer is f b , the static zero bias of the accelerometer and gyroscope is B a =(B ax ,B ay ,B az ) T and B ω =(B ωx ,B ωy ,B ωz ) T , the standard deviation of the accelerometer noise and the standard deviation of the gyroscope noise are and The sampling frequency is 1 / T s , then the error model of the inertial navigation source at time t is e INS ~N(μ INS ,Σ INS ).

[0142] in,

[0143]

[0144] Among them, [f n ×] indicates The antisymmetric matrix of the accelerometer measurements under the system, and Direction cosine matrix By G b OK, [f n ×].^2 indicates that [f n ×] squares each element of the matrix, and diag(·) is a diagonal matrix construction function.

[0145] Assume that after the inertial navigation correction, p n =(x,y,z) T At this moment, the positioning result of the drone to be positioned is p n =(x,y,z) T , then the information geometry probability model of the inertial navigation source is p INS ~N(p n -μ INS ,∑ INS ).

[0146] The drone is equipped with an ADS-B navigation system. Assume that the position of the node to be located is p = (x, y, z) T , the number of communicable ADS-B base stations is n (n ≥ 4), and the position of the i-th base station is P i ={(X i ,Y i ,Z i ) T |i=1,2,3,…,n}, the distance information between the node to be located and the communicable base station is d=(d1,d2,…,d n ) T .

[0147] According to the Euclidean distance calculation formula, we have the following equations:

[0148]

[0149] In order to reduce the influence of the inherent error of the UAV to be positioned, the subtraction of two adjacent equations can be obtained:

[0150]

[0151] The above formula can be rewritten into matrix form:

[0152] Ap=b………………………………(34)

[0153] in,

[0154]

[0155] Let H = (A T A) -1 A T , its least squares solution is:

[0156]

[0157] The distance measurement information between different drones is set to be completely independent, and the position of the node to be located is obtained according to the least squares method. The covariance of

[0158]

[0159] in,

[0160]

[0161] Where D(·) is the variance function.

[0162] Assume that the random variable d satisfies d~N(μ,σ 2 ), then

[0163] D(d 2 )=2σ 4 +4μ 2 σ 2 ………………………………(39)

[0164] The location information μ of the node to be located RNS It can be solved by the least square method. The covariance matrix of the location information of the node to be located ∑ RNS It can be obtained by the above formula. Then the information geometry probability model of ADS-B navigation source can be expressed as p RNS ~N(μ RNS ,∑ RNS ).

[0165] Step 3: Compare the positioning results of the low-orbit satellite for the UAV A to be positioned obtained in step 1 and the collaborative UAV A obtained in step 2 i The relative positioning results of the target A are fused to obtain the fused positioning results, and a UAV cluster information geometric fusion positioning model is constructed.

[0166] In this embodiment, steps 1 and 2 obtain the positioning results of the low-orbit satellite, inertial navigation and ADS-B navigation system. The specific process of fusion positioning is as follows:

[0167] In the drone cluster of this embodiment, it is assumed that there are j positioning results from different sources of drones, and the probability density of the i-th positioning result is p i (x|θ), i=1,...,j, the probability density of the fusion positioning result is If p i , i = 1, 2, ..., j belong to the same probability distribution family S = {p(x|θ)}, the distribution family S is regarded as a Riemann manifold and has a Fisher information matrix. The optimization criterion is used to select the optimal point on the distribution family S The theoretical criteria for information fusion are obtained:

[0168]

[0169] Here, the distance d is the Fisher information distance between p1 and p2. Under the assumption of Gaussian distribution, the normal distribution with mean μ and variance ∑ is expressed as N(μ, ∑), μ f ,∑ f Represents the mean and covariance matrix of the fusion result, and rewrites the above formula as follows:

[0170]

[0171] The probability density distribution family S belongs to a kind of exponential distribution family, which has the following form:

[0172] p(x|θ)=exp{C(x)+θ T F(x)-ψ(θ)}……………………(42)

[0173] Where C(x) is a polynomial of x, F(x) is a sufficient statistic for the natural parameter θ, ψ(θ) is the potential function of the distribution, and the potential function of the zero-mean complex Gaussian distribution N(0, ∑) is:

[0174] ψ(θ)=-log(|∑|)………………………………(43)

[0175] Let R1 and R2 be any two points on the manifold S. According to the above two formulas, the KL divergence distance between R1 and R2 on the manifold is:

[0176]

[0177] Among them, tr(·) represents the trace of the matrix. KL divergence can be used as a measure of the information difference between matrices. It has information monotonicity, that is, the divergence value increases as the information difference between the two matrices becomes larger, so the divergence is regarded as a distance on the manifold.

[0178] For N probability density distributions p(μ i ,∑ i ), i = 1, ..., N, and merge them into a probability density p f (μ f ,∑ f ), the problem of finding the optimal fusion strategy for drone clusters becomes finding the minimum value of p(μ i ,∑ i ):

[0179]

[0180] Among them, μ i and∑ i Represents the mean and covariance matrix of the probability density distribution to be fused.

[0181] During the fusion process, the covariance matrix ∑ i , i=1,...,N are integrated, so that the problem is transformed into finding the minimum value of ∑ f :

[0182]

[0183] Further, the fusion positioning result μ is obtained f :

[0184]

[0185] Among them, ω i is the weight factor, 0≤ω i ≤1 and satisfy the constraints:

[0186] ω1+ω2+…+ω N =1………………………………(48)

[0187] The covariance matrix provides a scalar measure of the uncertainty of the estimate, expressing the largest linearly independent subset of the linear constraints as:

[0188] tr(∑ i )ω i -tr(∑ i+1 )ω i+1 =0,i=1,2,...,N-1………………(49)

[0189] make Combining the above two formulas, we get:

[0190]

[0191] The system matrix is ​​expressed as Weight factor ω i Through the Laplace expansion of the first and nth columns, the Laplace expansion of the first and nth columns of θ are:

[0192] The weight factor ω is derived from the above formula i Parsing the expression:

[0193]

[0194] The positioning accuracy simulation test of the UAV cluster positioning in the embodiment is carried out by using the method of the present invention. Figure 3As shown in the figure, the mean square error (MSE) is used to calculate the positioning accuracy, and compared with the current mainstream fusion methods including factor graph method, Kalman filter, and artificial neural network fusion method. The simulation test scene size is 5000m×5000m×500m, the number of drones in the cluster is n=10, and the communication distance between drones is L=2000m. Among them, the number of Monte Carlo simulations is T n =1000, number of iterations t n =50, UAV ranging error variance Error variance of pseudo-range measurement of low-orbit satellite experimental satellite

[0195] from Figure 3 It can be seen that when different fusion positioning methods are used for testing without prior position information of the drone cluster, the positioning convergence speed of the method of the present invention is the fastest. After receiving the low-orbit satellite experimental satellite signal and the drone ranging information, the method of the present invention only needs 3 iterations to achieve the convergence of the positioning result, and the fusion positioning result after convergence fluctuates very little over time; when the Kalman filter, artificial neural network, and least squares fusion methods are used, the number of iterations required for the fusion positioning convergence are 10, 13, and 4 times, respectively. After the least squares fusion method converges, the positioning still has obvious fluctuations. The above test results show that the method of the present invention can quickly process pseudorange information and ranging and direction finding information, realize the fusion of various navigation source positioning information in a short time, have significant advantages in real-time positioning, and can well realize the real-time positioning of drone clusters in a fast-moving state.

[0196] The positioning stability of the UAV cluster positioning in the embodiment is simulated and tested in the face of mutation error using the method of the present invention. Figure 4 As shown in the figure, the error mutation simulation test conducts three mutation error tests: the first test is conducted at the 5th iteration after the positioning algorithm is run, and the host computer randomly shuts down the low-orbit satellite positioning data link of 2 UAVs; the second test is conducted at the 15th iteration, and the low-orbit satellite positioning data link of 3 UAVs is randomly shut down; the third test is conducted at the 30th iteration, and the low-orbit satellite positioning data link of 4 UAVs is randomly shut down.

[0197] Depend on Figure 4It can be seen that in the three mutation error tests, mutation errors occurred in the positioning results of the drones due to the loss of the low-orbit satellite experimental satellite signal by some drones. Among the four methods, the error increment of the method of the present invention is the smallest. In the three tests, the errors increased by 0.42m, 0.58m, and 0.55m respectively. After the mutation error occurs, the positioning error of the present invention can be stabilized again only after 3-4 moments and return to the normal value. Since the number of drones that lost the low-orbit satellite signal in the third mutation error test is the largest, the mutation error generated by positioning is also the largest. The mutation errors based on artificial neural network, Kalman filter, and least squares positioning are 1.33m, 1.18m, and 1.35m, respectively. By comparison, it can be seen that the mutation error of the method of the present invention is the smallest, and the present invention has the fastest convergence speed after the mutation error occurs. The test results show that the present invention converts the positioning information of low-orbit satellites, inertial navigation, and ADS-B navigation systems into information geometry density functions, and realizes rapid fusion through KL divergence, which can effectively suppress the mutation of drone positioning results in occluded scenes, and can meet the continuous and real-time high-precision positioning requirements of large-scale drone clusters in occluded environments.

[0198] The method of the present invention is used to perform simulation tests on drone clusters of different sizes in the embodiment. Figure 5 As shown in the figure, the larger the size of the drone cluster, the denser the distribution of drones in the three-dimensional space, and the larger the number of cooperative drones for each drone, the richer the drone cooperative positioning. Due to the fast movement speed of low-orbit satellites and the small satellite coverage area, the visible low-orbit satellites of drones are changing rapidly, and it is difficult to maintain a satellite communication link for a long time. Therefore, for the cooperative positioning of drone clusters, the positioning effect of low-orbit satellites becomes better as the number of drones in the cluster increases.

[0199] like Figure 5 As shown, during the test, the scale of the drone cluster is changed every 20 iterations, and the number of iterations in the time periods of 0-20, 20-40, 40-60, 60-80, and 80-100 corresponds to the number of drones participating in the fusion of 10, 8, 6, 8, and 10. As the number of drones in the cluster changes, the loss or increase of navigation signals will affect the positioning of the drone cluster. Among them, the least squares method has the greatest impact, with a maximum error of 5.8m. Due to the high computational complexity, the artificial neural network has a very slow convergence speed, and the positioning error after convergence can reach 1.4m. The impact on the present invention is the smallest. When the scale of the drone cluster changes, the convergence speed is the fastest, and the mutation error generated fluctuates within 0.2m. The test results show that the present invention can adapt to the network topology under different drone densities, and has high positioning accuracy and mutation suppression capabilities.

[0200] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention without departing from the principles and intent of the present invention.

Claims

1. A method for collaborative positioning of UAV clusters and low-orbit satellites based on information geometry, characterized in that: The following steps are involved: Step 1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory, ensure that the number of pseudo-ranges of drones in the cluster is 4, and solve the three-dimensional covariance matrix of the relative positioning error: Step 1-1: Construct a pseudo-range reuse algorithm model between drones based on factor graph theory: equip each drone in the cluster with a low-orbit satellite receiver, ranging equipment, inertial navigation, and ADS-B navigation to ensure that each navigation system can work independently. The drone can communicate with the low-orbit satellite and establish communication with other drones within the communication distance, and then measure the current distance and direction to the other party; Step 1-2: Process of obtaining the positioning of low-orbit satellites: Each drone in the drone cluster communicates with the low-orbit satellite through its own low-orbit satellite receiver, and obtains the low-orbit satellite navigation message and pseudo-range data through pseudo-range multiplexing of surrounding drones, constructs a positioning solution equation group based on the low-orbit satellite pseudo-range, and uses the least squares method and Newton iteration method to solve the equation group to solve the low-orbit satellite positioning coordinates of each drone in the cluster; Step 1-3: According to the positioning error introduced by relative positioning during pseudorange reuse, the three-dimensional covariance matrix of the low-orbit satellite positioning result of each UAV is obtained by analyzing the error transmission model of the nodes in the factor graph; Step 2: According to the navigation principles and error transfer models of low-orbit satellites, inertial navigation and ADS-B navigation, solve the positioning results and their three-dimensional covariance matrices of low-orbit satellites, inertial navigation and ADS-B navigation; It is assumed that UAV A1 has n visible low-orbit satellites, and the signal-to-noise ratio of each low-orbit satellite is SN i , i = 1,…,n, the low-orbit satellite signal with the maximum signal-to-noise ratio is defined as the optimal, then the optimal pseudorange selection criterion is: [min_SN,i]=max(SN i ) Among them, min_SN is the signal-to-noise ratio of the optimal low-orbit satellite, i is the serial number of the low-orbit satellite, and the optimal pseudorange of UAV A1 comes from the low-orbit satellite L i ; Assuming that the number of drones that drone A1 can communicate with is m, drone A1 can indirectly use the pseudo-range data of the communicative drones. There are m+1 pseudo-range data that drone A1 can use. All pseudo-range data are constructed into the pseudo-range data matrix of drone A1. Through the minimum GDOP criterion of pseudo-range, four matrices with the best geometric properties and the most suitable for drone A1 positioning are selected from the pseudo-range data matrix to solve the low-orbit satellite positioning results. Step 3: Compare the positioning results of the low-orbit satellite for the UAV A to be positioned obtained in step 1 and the collaborative UAV A obtained in step 2 i The relative positioning results of the target A are fused to obtain the fused positioning results, and a UAV cluster information geometric fusion positioning model is constructed.

2. The method for cooperative positioning of UAV clusters and low-orbit satellites based on information geometry as claimed in claim 1, characterized in that: In step 1, any target node m to be located q As an example, an internal factor graph is established. Since the coordinate axes are independent of each other in the coordinate solution process of the node to be located, the three-dimensional positioning problem is transformed into a one-dimensional problem on three branches, which is divided into x-coordinate, y-coordinate and z-coordinate. The three parts are connected by Euler's formula and expressed as: The constraint relationship formula is: in, Δx iq , Δy iq and Δz iq They represent the relative distances between the node to be located q and its neighbor node i on the coordinate axis respectively; d iq Indicates the relative distance between nodes; (x i ,y i ,z i ) represents the coordinate value of neighbor node i; (x q ,y q ,z q ) represents node m q The coordinates of .

3. The method for cooperative positioning of UAV clusters and low-orbit satellites based on information geometry as claimed in claim 1, characterized in that: In step 2, the positioning result of the low-orbit satellite is: in, x k , x k-1 are the calculation results of the kth and k-1th Newton iterations respectively; Δx, Δy, Δz, Δτ are the differences between the two solutions in the three directions during the Newton iteration process and the differences between the two clock error solutions. τ is the clock difference between the low-orbit satellite and the ground drone.

4. The method for cooperative positioning of UAV clusters and low-orbit satellites based on information geometry as claimed in claim 1, characterized in that: In step 2, it is assumed that at time t after the inertial navigation correction, the positioning result of the UAV to be positioned is p n =(x,y,z) T , then the information geometry probability model of the inertial navigation source is: p INS ~N(p n -m INS ,∑ INS ) in, in, B a =(B ax ,B ay ,B az ) T is the static zero bias of the accelerometer; B ω =(B ωx ,B ωy ,B ωz ) T is the static bias of the gyroscope; is the standard deviation of the accelerometer noise; is the standard deviation of the gyroscope noise; 1 / T s is the sampling frequency; [f n ×] is The antisymmetric matrix of the accelerometer measurements under the system, and Direction cosine matrix By G b Sure; [f n ×].^2 indicates that [f n ×] Square each element of the matrix; diag(·) is a diagonal matrix construction function.

5. The method for cooperative positioning of UAV clusters and low-orbit satellites based on information geometry as claimed in claim 1, characterized in that: In step 2, the information geometry probability model of the ADS-B navigation positioning source can be expressed as: p RNS ~N(μ RNS ,∑ RNS ) in, μ RNS is the location information of the node to be located; ∑ RNS is the covariance matrix of the location information of the node to be located.

6. The method for cooperative positioning of UAV clusters and low-orbit satellites based on information geometry as claimed in claim 1, characterized in that: In step 3, the process of obtaining the fusion positioning result is: The positioning results from low-orbit satellites, inertial navigation and ADS-B navigation are constructed as a three-dimensional Gaussian probability density function. Combined with the Riemann manifold under the Gaussian distribution family, the Fisher information distance is approximated by calculating the KL divergence of the probability density function of the positioning result on the Riemann manifold. When the sum of the Fisher information distances from a point in the Riemann manifold to each information source is the shortest, the probability density function corresponding to the point is the fused positioning result, thereby achieving precise positioning.

7. The method for cooperative positioning of UAV cluster low-orbit satellites based on information geometry according to claim 1 or 6, characterized in that: In step 3, the fusion positioning result is: in, ω i is the weight factor, 0≤ω i ≤1 and satisfy the constraints ω1+ω2+…+ω N =1; μ i and∑ i Represents the mean and covariance matrix of the probability density distribution to be fused; μ f ,∑ f Represents the mean and covariance matrix of the fusion result.

Images