Unmanned aerial vehicle target positioning method based on improved variation firefly optimization particle filtering

Abstract

The invention relates to an unmanned aerial vehicle target positioning method based on improved mutation firefly optimization particle filtering, and the method comprises the steps: fusing an optimization mechanism of a firefly algorithm in a frame of a particle filtering algorithm, and guiding particles to move towards a high-likelihood region; a multi-variation strategy and an elastic force mechanism are introduced to change a particle interaction mode, and the problem of particle degradation caused by serious nonlinear factors and excessive optimization is solved; the target position is optimally estimated by using the mutation firefly optimization particle filter, the number of particles required by the operation of the standard particle filter algorithm is reduced, and the robustness and the positioning precision of the algorithm are improved.

Description

Technical Field

[0001] The present invention relates to the field of unmanned aerial vehicle (UAV) target positioning methods, and in particular to an unmanned aerial vehicle (UAV) target positioning method based on improved variant firefly optimized particle filtering. Background Art

[0002] The application of optoelectronic platforms on UAVs is one of the core elements of their intelligent development. With the deep integration of computer vision and optoelectronic sensing technology, optoelectronic platforms carried by UAVs have made significant technological breakthroughs in target recognition, regional search, and remote sensing mapping, making UAVs play a vital role in key areas such as transportation, aviation, and electricity.

[0003] Using optoelectronic platforms to accurately locate and estimate the state of targets is a prerequisite for UAVs to perform missions. However, there are many factors that affect the target positioning of UAVs, and the transmission chain is long. In addition, in the actual mission execution process, they are often affected by harsh environments such as strong winds and vibrations, which causes the airborne optoelectronic platform to be seriously affected by nonlinear factors during mission observation. In addition, when UAVs perform reconnaissance and surveillance missions in dangerous target areas or non-airspace areas, optoelectronic platforms mainly use high-altitude oblique vision to image large areas at long distances. The observation conditions often have extreme characteristics of large inclination angles, long distances, long focal lengths, and strong nonlinearity, which causes the observation equations of optoelectronic platforms to be seriously ill-conditioned, and the target positioning accuracy is very sensitive to the observation errors caused by disturbances. Therefore, it is of great significance to carry out research on target positioning methods with strong computational robustness and high positioning accuracy.

[0004] In order to solve the problem of significantly reduced target positioning accuracy due to the above-mentioned strong nonlinear factors, scholars at home and abroad have carried out a series of work mainly from the two paths of image matching positioning and geometric positioning. Among them, image matching positioning is to register the measured image and the reference image based on the image affine invariance when the reference image of the target area is known, so as to eliminate the modal variance and geometric deformation inherent in image matching. In 2004, David proposed a scale-invariant feature transform (SIFT) method, which is mainly used to perform reliable registration between different perspective images of objects or scenes, and has become one of the most popular solutions for image matching. The FSG-SIFT algorithm proposed by Change et al. introduced feature slope calculation and feature point grouping on the basis of the scale-invariant feature transform (SIFT) algorithm, and used random sample consensus to refine the matching, which effectively improved the feature matching accuracy. Yang et al. proposed a method based on feature registration and phase-correlation registration, which retains the high versatility of feature registration and the high accuracy of extended phase-correlation registration, and has high robustness and accuracy. Based on the SIFT method, Paul et al. used the scale, direction and translation differences between the input images to remove outliers, thereby retaining more correct matching pairs. Mohammadi et al. used a method based on self-similarity measurement to compress robust image features into feature descriptors, resist the difference in illumination between images, and improve the accuracy and reliability of multimodal remote sensing image matching under significant time differences. Yao et al. improved the traditional image feature matching method and used a co-occurrence filter to construct a new co-occurrence scale space, thereby extracting feature points in the new scale space and effectively improving nonlinear distortion. Tian Dan et al. combined the ORB feature extraction operator with the BRISK feature description operator based on the idea of ​​consistent local structure of similar images. The proposed method effectively solved the problems of high mismatch rate and slow matching speed in image matching.

[0005] However, under extreme observation conditions, the resolution and clarity of aerial images are greatly reduced, and they show obvious perspective transformation effects. The image matching positioning method using the geographic location information of the reference map has low accuracy and is difficult to use for real-time positioning and positioning of moving targets. Geometric positioning uses optoelectronic platform and UAV sensor data as observation quantities, constructs observation equations and state equations based on the principle of photogrammetry, does not require prior knowledge and has good real-time performance, and has become the mainstream of UAV target positioning methods. According to whether the optoelectronic platform actively emits radiation to the target during observation, geometric positioning methods can be mainly divided into active positioning and passive positioning. Active positioning uses a laser target indicator to obtain the distance information between the UAV and the target, and then locates the target in combination with the angle information, with high positioning accuracy. Chen Danqi et al. proposed a Gauss-Newton iterative nonlinear target positioning algorithm based on laser ranging values, which can meet the requirements of positioning accuracy and real-time performance at the same time. However, the active positioning method actively emits radiation outward, which usually exposes the position of the UAV itself during the mission. In addition, since the laser ray coincides with the platform optical axis, active positioning can only obtain the distance of the target at the center of the field of view, and cannot locate multiple targets at the same time. Passive positioning methods passively receive the signal source radiated by the target, and can make the best estimate of the target's geographic coordinates through the corresponding algorithm. They are widely used in the field of UAV positioning: Bai Guanbing et al. used Kalman filter (KF) based on the geometric mathematical model to achieve real-time positioning of multiple targets in the same field of view. However, the highly nonlinear characteristics of the optoelectronic system observation model make it difficult for the traditional KF algorithm to meet the actual accuracy requirements. Although KF derivative algorithms such as Extended Kalman Filter (EKF), Cubic Kalman Filter (CKF), and Unscented Kalman Filter (UKF) can handle nonlinear problems to a certain extent, the above algorithms are all based on Gaussian noise, and their ability to solve positioning problems with strong nonlinear characteristics is relatively limited. The particle filter (PF) algorithm does not require linear assumptions on the system model and observation model, and has better universality. However, the PF algorithm also has its shortcomings. As the algorithm continues to run, the accumulated nonlinear factors will cause the loss of sample validity and diversity, and will inevitably produce particle degradation, seriously affecting the accuracy of passive positioning of UAV targets. Summary of the invention

[0006] In view of the serious nonlinear factors affecting the positioning of UAV targets, this paper proposes a UAV target positioning method based on Improved Mutant Firefly Algorithm (IMFA-PF). The main idea of ​​this method is to use the state equation of the imaging system of the UAV optoelectronic platform to optimally estimate the geographical location of the target by combining swarm intelligence optimization and particle filtering technology, thereby achieving accurate positioning of ground targets by UAV under serious nonlinear influences.

[0007] The present invention adopts the following technical solutions:

[0008] The UAV target positioning method based on improved mutant firefly optimized particle filter includes the following steps:

[0009] Step 1: Collect the real-time position observation value of the target, add a random noise to the real-time observation value each time, repeat N times, and obtain N position data. Each position data corresponds to the state of a particle. N particles form a particle swarm. The position state corresponding to the particle swarm at time k-1 constitutes the particle state set

[0010] Step 2: According to the particle state set at time k-1 Using the state equation of the positioning system to predict the particle state set at time k Calculate the luminous brightness of each particle and select the particle with the largest brightness as the optimal particle. The state value corresponding to the optimal particle is

[0011] Step 3: Calculate the optimal particle Attraction to other particles The position of each particle at time k is updated and calculated according to the attraction degree. The updated position of the i-th particle is Repeat the position swap mutation for m1 times on the updated particle swarm, retain the particle with the best position each time, and obtain the swapped mutated particles. The state set corresponding to the swapped mutated particles is

[0012] Step 4: Repeat the position differential mutation for m2 times on the particle swarm after the exchange mutation to form a new particle swarm after the differential mutation. Select the N particles with the best performance from the particle swarm before the position differential mutation and the new particle swarm after the differential mutation as the particle swarm for the next iteration. The state set corresponding to the particle swarm is

[0013] Step 5: Calculate the density ρ of the particle group obtained in the previous step k , judge ρ k Is it greater than the preset threshold ρ th, then the particle group is dispersed according to the elastic mechanism to obtain the dispersed particle group. The position state of the i-th dispersed particle is Otherwise, keep the position of the particle group unchanged;

[0014] Step 6: Use the target observation value z at the current k moment k Calculate normalized particle weights

[0015] Step 7: Output the target location estimate at time k

[0016]

[0017] Furthermore, the formula for updating the position of each particle at time k according to the attraction in step 3 is:

[0018]

[0019] In the formula, is the updated state value of particle i after moving at time k; is the initial state value of particle i at time k, α∈[0,1] is the step factor, and rand is a random number uniformly distributed on [0,1];

[0020] The calculation formula of α is

[0021]

[0022] Where c is the scale factor and μ is the position factor.

[0023] Furthermore, the specific method of position swap mutation in step 3 is:

[0024] A particle is randomly selected and the two dimensions of its state space are exchanged to generate a new particle. The particle with the best position is retained in the new particle and the original particle. One position exchange mutation is completed, and the mutation operation is repeated m1 times.

[0025] Furthermore, the specific method of each position difference mutation in step 4 is:

[0026] Randomly select two particles to perform differential mutation to obtain new particles. and Perform position differential mutation to obtain new particles The calculation formula is:

[0027]

[0028] Where F is the scaling factor.

[0029] Furthermore, in step 5, ρ k The calculation formula is:

[0030]

[0031] Where num is the operation for calculating the number of particles in the judgment area; r j is the radius of the judgment area; N is the total number of particles;

[0032] The calculation formula is:

[0033]

[0034] In the formula, ξ is the rebound coefficient, k is the elastic coefficient, For elasticity, is the compression amount of particle i at time k.

[0035] After adopting the above technical solution, the present invention has the following advantages compared with the prior art:

[0036] (1) The optimization mechanism of the firefly algorithm is integrated into the framework of the particle filter algorithm to guide particles to move to the high likelihood area;

[0037] (2) Introducing multiple mutation strategies and elastic force mechanisms to change the mode of particle interaction and solve the problem of particle degradation caused by severe nonlinear factors and over-optimization;

[0038] (3) Use the mutant firefly optimized particle filter to optimally estimate the target position, reduce the number of particles required to run the standard particle filter algorithm, and improve the robustness and positioning accuracy of the algorithm.

[0039] The present invention is described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 It is a schematic diagram of the basic coordinate system;

[0041] Figure 2 It is the curve diagram of particle attraction variation with different light intensity attenuation coefficients;

[0042] Figure 3 This is a comparison diagram of the system filter state estimation;

[0043] Figure 4 This is a comparison chart of the absolute values ​​of the estimation errors;

[0044] Figure 5 Particle distribution diagram at different times;

[0045] Figure 6 This is a comparison chart of the measured results of PF and IMFA-PF algorithms (N=50);

[0046] Figure 7 Comparison of measured results between PF and IMFA-PF algorithms (N=500). DETAILED DESCRIPTION

[0047] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0048] 1. Principle of UAV target positioning

[0049] 1.1 Basic coordinate system

[0050] The process of UAV target positioning mainly involves 7 basic coordinate systems, such as Figure 1 shown.

[0051] 1) Geodetic coordinate system GCF(O gcf -LBH)

[0052] Figure 1 The GCF coordinate system in (a) takes the center of the Earth ellipsoid as its origin and is a three-dimensional curvilinear coordinate system that uses longitude L, latitude B, and altitude H to describe the position of points on the Earth's surface.

[0053] 2) Geocentric coordinate system ECEF(O ecef -X ecef Y ecef Z ecef )

[0054] Figure 1 The ECEF coordinate system in (a) also takes the center of the Earth ellipsoid as its origin, and X ecef The positive direction of the Y axis points to the North Pole of the Earth. ecef The positive direction of the axis points to the intersection of the prime meridian and the equator, Z ecef The positive direction of the axis is perpendicular to the other two axes and points toward the equator.

[0055] 3) Navigation coordinate system NED(O ned -X ned Y ned Z ned )

[0056] Figure 1 (a) The black coordinate system with its origin at point K and Figure 1 (b) The green coordinate system with the origin at the drone is the same coordinate system, namely the NED coordinate system. NED uses the center of mass of the drone as the origin of the reference system, and X ned The positive direction of the axis points to the north of the earth, Y ned The positive direction points to the east of the earth, Z ned The positive direction of the axis is perpendicular to the earth's surface and points downward.

[0057] 4) Body coordinate system HRD(O hrd -X hrd Y hrd Z hrd )

[0058] Figure 1 (b) The black coordinate system with the origin at the drone is the HRD coordinate system, whose origin coincides with NED. hrd The positive direction of the axis points to the front of the aircraft, hrd The positive direction of the axis points to the right wing of the aircraft, Z hrd The positive direction of the axis is downward along the vertical axis of the carrier.

[0059] 5) Camera coordinate system C(O) c -X c Y c Z c )

[0060] Figure 1 (b) The orange coordinate system with its origin at the pod is the C coordinate system, whose origin is the optical center of the camera, X c Axis positive direction and Y c The positive direction of the axis is parallel to the image plane, pointing to the right side of the camera and the lower side. c The positive direction of the axis points to the imaging direction along the optical axis.

[0061] 6) Image physical coordinate system W(ox w y w )

[0062] Figure 1 (b) The black coordinate system with its origin at the center of the image plane is the W coordinate system, x w Axis and y w The X axis c Axis and Y c Axis parallel.

[0063] 7) Image coordinate system I (I-xy)

[0064] Figure 1 (b) The black coordinate system with its origin at the upper left corner of the image plane is the I coordinate system, and the x-axis and y-axis are also the same as the X-axis. c Axis and Y c Axis parallel.

[0065] 1.2 State equation and observation equation of UAV positioning system

[0066] 1) State equation of target positioning The purpose of target positioning is to estimate the position information of the observed object. Therefore, the position of the target in the navigation coordinate system is selected as the state vector of the system, that is,

[0067] x k =[xk ,y k ,z k ] T (1)

[0068] Where: x k ,y k and z k are the coordinates of the target in the navigation coordinate system at time k. Then the state transition model of the system can be expressed as

[0069] x k =f(x k-1 ,w k-1 )=Φ k|k-1 x k-1 +w k-1 (2)

[0070] Where: subscripts k and k-1 represent time k and time k-1 respectively; Φ k|k-1 represents the state transition matrix from time k-1 to time k. For a stationary target, it is a 3×3 unit matrix. k-1 is the system process noise.

[0071] 2) Observation equation for target positioning

[0072] The observation equation of the positioning system can transform the pixel coordinates in the I coordinate system to the earth coordinate system through coordinate transformation. Assume that the system observation at time k is the actual pixel coordinate z of the target in the I coordinate system k =[u,v] T , which can be obtained by computer vision tracker. Then the observation equation is

[0073] z k =g(x k ,v k )=h(x k )+v k (3)

[0074] Where: v k is the observation noise; h(·) is the nonlinear function of the optoelectronic platform to convert the target position into pixel coordinates. The specific solution process is as follows:

[0075] The coordinates of the target P in the C coordinate system are (x c ,y c ,z c ), from the collinear equation of photogrammetry:

[0076]

[0077] In the formula, u′ k and v′ kare the observed quantities calculated by the state vector; m1~m 12 The elements of the projection matrix are determined by the intrinsic and extrinsic parameters of the camera. After obtaining the coordinates of the target in the C coordinate system, they are converted to the geodetic coordinate system through the earth ellipsoid model.

[0078] 2 Based on the improved firefly optimized particle filter method

[0079] 2.1 Particle Filter Principle

[0080] Particle filtering is a method that uses the theoretical basis of Monte Carlo simulation and the law of large numbers to replace the traditional integral operation with the sample mean to effectively obtain the minimum variance estimation method of the system state. Its main concept is to use N independent and identically distributed particles and their corresponding weights To approximate the posterior probability distribution function. Assuming that the initial posterior probability density of the positioning system state is p(x0|z0)=p(x0), the state prediction equation is

[0081] p(x k |z 1:k-1 )=∫p(x k |x k-1 )p(x k-1 |z 1:k-1 )dx k-1 (5)

[0082] The state update equation is

[0083]

[0084] The posterior probability density function is estimated as

[0085]

[0086] Where δ(·) is the Dirac function, is the importance weight of the i-th particle at time k, and its formula is

[0087]

[0088] Finally, normalize the weights and output the final filter state

[0089]

[0090] 2.2 Firefly Algorithm

[0091] The Firefly Algorithm (FA) is a heuristic swarm optimization algorithm, which is inspired by the behavior of fireflies in nature. In nature, fireflies attract mates or prey through their luminous properties, and their brightness is usually proportional to the individual's luminous brightness and attractiveness. The specific description of the FA algorithm is as follows:

[0092] 1) Relative brightness of fireflies I

[0093]

[0094] Where I0 is the maximum brightness of the firefly; γ is the light intensity attenuation coefficient, which controls the speed at which the light intensity decays with distance; r ij is the Euclidean distance between firefly i and firefly j.

[0095] 2) Firefly attraction β

[0096]

[0097] Where β0 is the maximum attraction, which is usually set to 1.

[0098] 3) The position update formula for each firefly moving towards other brighter fireflies is

[0099] x i =x i +β×(x j -x i )+α×( rand -0.5) (13)

[0100] In the formula, x i 、x j is the spatial position of fireflies i and j; α∈[0,1] is the step size factor, and rand is a random number uniformly distributed on [0,1] to prevent the algorithm from falling into the local optimum too early.

[0101] 2.3 Firefly Optimized Particle Filter Based on Multiple Mutations

[0102] Directly fusing FA with particle filtering will easily lead to local optimality and particle degradation. In order to improve particle diversity and algorithm convergence accuracy, a firefly optimized particle filter based on multiple mutation mechanism is proposed. First, the initialized particles are regarded as fireflies, and the distribution of particles is optimized by applying attraction and movement mechanisms, and the particles are guided to migrate to high likelihood areas; then, multiple mutation strategies are introduced to change the mode of particle interaction and explore new solution spaces to solve the nonlinear problem in the observation process; finally, the distribution density of particles in space after each iteration is detected. If the density is too high, an elastic mechanism is established between the optimal particle and the particles in the surrounding area to disperse the over-concentrated particles in the area, thereby achieving a more reasonable particle distribution and preventing over-optimization of the algorithm.

[0103] 2.3.1 Particle attraction and movement mechanism

[0104] 1) Particle luminescence brightness

[0105] Using the latest measurement value z of the filter at time k k Calculate the luminous brightness of particle i To evaluate the quality of the particle position, the specific formula is

[0106]

[0107] In the formula, is the predicted value of particle i at time k; ||·||2 is the Euclidean distance between the measured value and the predicted value of particle i.

[0108] The particle with the highest luminous brightness has the most advantageous position, so the optimal particle is the particle with the highest luminous brightness at time k. The specific calculation formula is:

[0109]

[0110] In the formula, is the state value of the ith particle at time k; I(·) is the operation to calculate the luminous brightness.

[0111] 2) Particle attraction

[0112] Particle i moves towards the optimal particle P best The distance moved depends directly on the attraction between the two, and the formula is

[0113]

[0114] In the formula, for k The attraction of the optimal particle to particle i at the moment.

[0115] 3) Particle movement

[0116] The attraction of the optimal particle to particle i can guide its state movement, and the position update formula after movement is:

[0117]

[0118] In the formula, is the updated state value of particle i after moving at time k; is the initial state value of particle i at time k.

[0119] 2.3.2 Multiple mutation mechanisms

[0120] In order to solve the problem of unreasonable particle aggregation and low-quality convergence caused by over-optimization, a multiple mutation mechanism is introduced to change the interaction between particles.

[0121] 1) Variation of step size α

[0122] In order to avoid particles oscillating near the optimal solution in the later stage of the algorithm, Levy variation is introduced in the step size α, that is, the non-negative Levy continuous probability distribution is used to replace the original randomization parameter. k The calculation formula of the time step α is:

[0123]

[0124] Where c is the scaling factor, which mainly affects the spread of the distribution and takes a value of 0.5; μ is the position factor, which mainly affects the left-right translation of the distribution and takes a value of 1.2. In the early stages of the algorithm, a larger step size α can expand the search range of the particles, enhance their global search capabilities, and prevent the algorithm from converging prematurely. As the iteration process proceeds, the step size gradually decreases, so that in the later stages of the algorithm, a smaller step size helps the particles to explore more deeply in the state space, thereby more effectively seeking the optimal solution.

[0125] 2) Position swap mutation

[0126] After any particle moves to the latest position, two dimensions of the state space are randomly selected for exchange to generate a new particle. Then the particle with the best position is always retained, and the mutation operation is repeated m1 times.

[0127] 3) Position differential variation

[0128] Randomly select two particles at time k and Perform position differential mutation, the calculation formula is

[0129]

[0130] Where F is the scaling factor, which is 0.5. The position differential mutation operation is performed m2 times to generate a new population. Then, the N particles with the best performance are selected from the original population and the newly generated population as the population for the next iteration. The purpose of position differential mutation is to explore new solution spaces by introducing mutation operations, thereby helping the algorithm to better search for the global optimal solution.

[0131] 2.3.3 Elastic force mechanism of high-density particles

[0132] After the particles complete an iteration and mutation, the density of the remaining particles around the optimal particle needs to be detected to achieve a more reasonable distribution of particles in the state space, ensuring that most particles are concentrated around the optimal particle, while ensuring that other areas also have appropriate particle distribution.

[0133] 1) Particle density determination

[0134] The particle density ρ at time k k for

[0135]

[0136] Where num is the operation for calculating the number of particles in the judgment area; r j is the radius of the judgment area; N is the total number of particles.

[0137] 2) Particle position update

[0138] If the particle density ρ k Exceeding the preset threshold value th , indicating that the distribution of particles around the optimal particle is too dense and needs to be adjusted. Therefore, an elastic mechanism is established between the optimal particle and the surrounding particles. This mechanism ensures that the distance between particles is proportional to their compression degree, that is, the greater the compression, the stronger the elastic force generated. Then the compression of particle i at time k is for

[0139]

[0140] According to Hooke's theorem, the corresponding elastic force is

[0141]

[0142] Where k is the elastic coefficient. The distance of the rebound should be proportional to the degree of compression, but this distance needs to be kept below the particle compression. So the rebound coefficient is defined as

[0143]

[0144] The rebound coefficient ξ drives the particle to bounce in the opposite direction of compression, and the position update formula of the bounced particle i is:

[0145]

[0146] 2.4 Specific implementation steps of the algorithm

[0147] The main idea of ​​this algorithm is to optimize the particle distribution through the attraction and movement mechanism of fireflies. At the same time, multiple mutation strategies and elastic force mechanisms are introduced to further optimize the particle interaction mode and solve the defects of the algorithm falling into local optimality and particle degradation. The specific process of the algorithm is as follows:

[0148] Step 1: When k=0, set the relevant important parameters: maximum attraction β0, light intensity attenuation coefficient γ, scale coefficient c, position coefficient μ, zoom coefficient F, judgment area radius r j , density threshold ρ th , elastic coefficient k, number of particle swarms N, initialize sample set And obey the importance density function:

[0149]

[0150] Step 2: Predict the particle state at time k The particle set at time k-1 And the state equation of the positioning system (Equation (2) and Equation (4)), we can predict

[0151] Step 3: Set the weight of each particle to w k =1 / N.

[0152] Step 4: Calculate the brightness of the particles And update the global optimal particle The improved firefly search algorithm is used to evolve particles. The luminous brightness of each particle is calculated according to formula (14), and the optimal particle is updated by formula (15).

[0153] Step 5: Calculate particle attraction And update the particle position: Calculate the optimal particle by formula (16) Attraction to other particles Then, the position of the particle after the move is obtained by equation (17). After each particle moves, the position swap mutation is repeated m1 times, and the relatively brightest firefly is saved. The updated particle position is

[0154] Step 6: Update the particle swarm Perform position differential mutation: Perform position differential mutation on the particle swarm according to formula (19) to explore new solution space and form a new particle swarm after mutation.

[0155] Step 7: Particle density detection: Calculate the particle density ρ within the radius using formula (20) k , if it is greater than the preset threshold ρ th , according to equations (21) to (24), disperse the particle swarm The updated particle swarm position is

[0156] Step 8: Normalized weight calculation: According to equations (9) and (10), using the latest observation value z k Calculate normalized particle weights

[0157] Step 9: Output target location estimate:

[0158]

[0159] 2.5 Algorithm time complexity analysis

[0160] IMFA-PF replaces the system resampling process of the traditional particle filter through the iterative optimization process of the heuristic algorithm. In order to evaluate the time cost brought by the algorithm improvement, the time complexity of the resampling process of the PF algorithm and the optimization process of IMFA-PF are analyzed respectively. Assume that the number of particles in both algorithms is N and the number of algorithm iterations is T.

[0161] 1) Time complexity analysis of PF resampling process

[0162] In the resampling phase of PF, the system resampling process of each iteration mainly involves generating a random starting point and traversing the weights at fixed intervals, so its complexity is O(N). The total time complexity of the PF algorithm initialization and resampling process is O(N×T).

[0163] 2) IMFA-PF time complexity analysis

[0164] In the iterative optimization stage of IMFA-PF, the state update of the particle swarm mainly involves three steps: calculation of luminous brightness, attractiveness, and position update. Each particle involves m1 position swap mutations when updating its position, so the time complexity is O((3+m1)N). After the update, m2 position differential mutation operations are performed, with a time complexity of O(m2N). After that, density detection is introduced to adjust the state of the high-density particle swarm, with a time complexity of O(2N). Since high-density detection is performed only when the algorithm is over-optimized, the number of executions should be less than the number of iterations T. Therefore, the total time complexity of the IMFA-PF algorithm initialization and resampling process should be less than O((5+m1+m2)×N×T).

[0165] In summary, compared with the traditional PF algorithm, the time complexity of the proposed algorithm has increased under the same conditions, but for the UAV observation system with a high degree of nonlinearity, the PF algorithm is very prone to particle degradation. The increased complexity of the proposed algorithm is mainly concentrated in the mutation operation, which effectively expands the search space of the optimal solution. If the standard PF algorithm is to achieve the same accuracy, more time and more particles are used, which is consistent with the conclusions of the simulation experiment.

[0166] 3 Experiments and results analysis

[0167] 3.1 Simulation Experiments and Results

[0168] In order to evaluate the specific performance of the IMFA-PF algorithm proposed in this invention, a highly nonlinear unstable system model is selected for simulation experiments. Its state equation and observation equation are:

[0169]

[0170] In the formula, e k ~N(0,Q) and v k ~N(0,R) are system noise and measurement noise respectively. The experimental platform is AMD R9-9500HX, 3.3GHz and 32GB RAM, and the simulation environment is Python 3.8.

[0171] 3.1.1 Important parameter settings

[0172] The important parameters in the filtering algorithm have a significant impact on its performance. Taking the light intensity attenuation coefficient γ as an example, we conduct an experimental study and analysis to obtain the optimal parameters. The light intensity attenuation coefficient determines the optimal particle at time k. For particle x kThe attraction directly affects the degree and efficiency of optimization. In order to test the influence of different light intensity attenuation coefficients on the attraction and obtain the optimal parameters that take into account both accuracy and efficiency, the light intensity attenuation coefficient γ is taken as 0.01, 0.05, 0.1, 0.5, 1, 1.5, and 2 respectively. The curve of attraction β changing with the distance r between the particle and the optimal particle when different γ values ​​are taken is as follows: Figure 2 shown.

[0173] Depend on Figure 2 It can be seen that at the same distance, as the light intensity attenuation coefficient increases, the optimal particle For particle x k The attraction gradually decays, and the larger the light intensity attenuation coefficient, the faster the attraction decays. When γ = 0.5, for particles with a distance greater than 3, the attraction of the optimal particle approaches 0. After one iteration, almost all particles with a radius less than 3 gather towards the optimal particle. This causes the particles to converge too quickly, which is not conducive to iterative optimization. When the value of γ increases, the decay rate of the attraction also accelerates, resulting in a faster convergence of particles. When γ = 2, after one iteration, all particles with a distance less than 1 will approach the optimal particle, which causes the algorithm to converge too quickly. When γ = 0.01, the optimal particle still has a certain attraction to particles with a distance greater than 10, and the optimization efficiency is low at this time. And for particles close to the optimal particle, if the attraction is too small, the effect of the mutation disturbance may be much greater than the attraction of the optimal particle, causing most particles to be in a disordered random search state, thereby affecting the fine search ability of the particles. When the light intensity attenuation coefficient is 0.05 or 0.1, as the distance between the optimal particle and other particles increases, the attraction increases relatively slowly, which is conducive to the gradual optimization of particles, and the particle optimization efficiency when γ = 0.1 is higher than that when γ = 0.05. Therefore, the light intensity attenuation coefficient γ is 0.1.

[0174] The other parameters of the IMFA-PF algorithm are tested in the same way. In addition, in order to evaluate the performance of the proposed algorithm for state estimation, particle filter, particle filter based on bat optimization algorithm (BA-PF), particle filter based on elastic firefly optimization algorithm (SFA-PF) and IMFA-PF are selected for state estimation performance comparison. The key parameters of the BAPF algorithm include attenuation coefficient α, enhancement coefficient λ, minimum pulse frequency f min and the maximum pulse frequency f max ; The key parameters involved in the SFA-PF algorithm include the maximum high likelihood particle ratio ε s , the maximum number of iterations is N i The specific settings of each algorithm parameter are shown in Table 1.

[0175] Table 1 Parameter settings of different algorithms

[0176]

[0177] 3.1.2 Algorithm performance test

[0178] 1) Particle noise resistance test

[0179] The accuracy of the filtering algorithm is related to the size of the noise. When the measurement noise is much larger than the system noise, the particle filter algorithm will face serious particle depletion problems. In the simulation experiment, the variance R of the measurement noise is set to 0.1, the variance Q of the state noise is set to 10 and 1 respectively, the total number of particle swarms is set to N = 100, and the number of filtering time steps is set to T = 100.

[0180] Figure 3 and Figure 4 The system state estimation results and estimated absolute errors of the four filtering methods under different state noise conditions are shown respectively. Comparing the two figures, it can be seen that when the variance of the state noise Q=10, the state noise is larger than the measurement noise, and the estimation error of PF is significantly larger than that of the other three algorithms. When Q=1, the state noise is lower, and the number of high errors and high risks is significantly reduced. This is because the larger the state noise, the farther the particle set is scattered, and the smaller measurement noise causes the high likelihood interval to narrow. In the early stage of the algorithm, the probability of particles with large state noise falling into this interval decreases. As the number of iterations increases, the particle weights in the high likelihood area decrease, resulting in most particles failing to fall into this interval many times, so that the system state estimation results gradually diverge and fall into a deeply degraded state. However, after optimization using the heuristic algorithm, BA-PF, SFA-FA and IMFA-PF can intervene in the state of the degraded particle group when the estimation results diverge, avoiding the transmission of erroneous particles to the next generation, thereby having a better tracking effect on the changes in the real state of the system.

[0181] 2) Particle diversity test

[0182] In order to verify the particle diversity during IMFA-PF filtering estimation, the variance of state noise and measurement noise is set to Q = 10 and R = 0.1 respectively, and the particle distribution of PF and IMFA-PF is compared, with a total number of particles N = 100. The particle distribution at time steps k = 10, k = 52, and k = 86 with large estimation errors are selected, as shown in Figure 2. Figure 5As shown in the figure. The black straight line is the real state of the system at the current moment, and the blue dots and red stars are the states of each particle of the PF and IMFA-PF algorithms respectively. As can be seen from the figure, at three different moments, the PF algorithm degrades due to the lack of particle diversity. The mean of the particle group deviates far from the real state, and the distribution is relatively scattered, which ultimately makes it difficult to ensure the filtering accuracy. In contrast, the IMFA-PF algorithm replaces the resampling process with iterative optimization of particles, and uses multiple mutation strategies and elastic mechanisms to expand the new solution space, thereby improving the particle distribution state. Therefore, most of the particles are concentrated near the real state, while a small number of particles deviate from the real value, which not only ensures the diversity of particles, but also has a higher precision state estimation at the current moment.

[0183] 3) Algorithm accuracy and real-time testing

[0184] In order to further study the accuracy and computational burden of the IMFA-PF algorithm, the number of particles N is set to 50, 100, and 200. The root mean square error (RMSE) is used to test the simulation system, and its calculation formula is:

[0185]

[0186] In the formula, x k is the true value of the state at time k; is the estimated value of the state at time k. The RMSE average, standard deviation and single-step average operation time of 50 simulations of the four algorithms are compared, and the results are shown in Table 2.

[0187] As shown in Table 2, with the increase in the number of particles involved in state estimation, the estimation accuracy of the four algorithms has improved to varying degrees, but the operating cost will also gradually increase. When the number of particles is 50, the error of the PF algorithm is larger than that of other algorithms, while when the number of particles is 200, the error gap between PF and other algorithms decreases, indicating that the accuracy of PF is more dependent on the number of particles. In addition, compared with PF, the RMSEmean and RMSEvar of the estimation results of the PF algorithm optimized by swarm intelligence are significantly reduced, which proves the effectiveness of swarm intelligence optimization in improving particle distribution and avoiding particle degradation. Under the condition of a small number of particles, IMFA-PF can better improve the particle degradation problem with the help of multiple mutation operations, so it can almost reach the accuracy of the other algorithms when the number of particles is large. Under the conditions of different numbers of particles, the attraction and movement mechanism of IMFA-PF ensures that the particle distribution is close to the high likelihood area, so that the RMSEvar is minimized, that is, the stability of state estimation is better than that of the other three algorithms. Secondly, the average single-step estimation time T of IMFA-PF and other algorithms is meanAt the same level, other algorithms need to further increase the number of particles to achieve the same accuracy and stability, which in turn increases the computational cost. Therefore, the IMFA-PF algorithm has better state estimation accuracy, and the number of particles and computational cost required to achieve a certain accuracy are less, meeting the requirements of filtering speed.

[0188] Table 2 Comparison of simulation results of different particle filter algorithms

[0189]

[0190]

[0191] 3.2 Actual flight experiment

[0192] In the UAV target positioning flight experiment, the angle measurement error of the UAV's roll angle, yaw angle, and pitch angle is less than 1.5 degrees; the angle measurement error of the optoelectronic pod's azimuth and elevation angle is less than 1 degree; the GPS positioning error of the UAV is less than 20cm. Before the UAV takes off, the ground station designs the flight route and plans to observe the target area. The UAV flies and observes autonomously according to the designed route, and the flight altitude is 600m. Figure 6 (c) As shown. During the observation of the ground target area, the measuring points are selected at equal time intervals. The GPS position data, IMU attitude data, pod angle measurement information and image data of each measuring point of the drone are sent to the ground station. Five 2×2m targets T1~T5 are set in the target area as measuring points.

[0193] In order to analyze the impact of the number of particles on the positioning of a single target T1, the total number of particle swarms N is set to 50 and 500 respectively, and the longitude and latitude estimation changes and positioning errors of the traditional PF and IMFA-PF algorithms are compared and analyzed. The variance of the measurement noise is set to R = 1.2, the variance of the state noise is set to Q = 0.0001, and the filtering time step is set to T = 100. The experimental results are shown in the figure. Figure 6 and Figure 7 The positioning error uses the circular probability error that can reflect the positioning accuracy of the target geographic location. The calculation formula is:

[0194]

[0195] In the formula, σ L and σ B are the positioning errors in longitude and latitude, R NT and R MT are the curvature radii of the meridian circle and the meridian circle, H T is the target height, B T is the target latitude.

[0196] Depend on Figure 6It can be seen that when the number of particles is 50, the positioning error of the PF algorithm decreases in the initial stage of filtering, but due to the accumulation of cumulative errors, the particles quickly degrade at time point 10, the positioning error gradually increases, and divergence occurs in subsequent estimates. However, IMFA-PF can use the swarm intelligence optimization algorithm to timely improve the impact of particle degradation, and further expand the solution space of particles with the help of mutation operations, and still maintain a good positioning effect under the condition of a small number of particles. Figure 7 It can be seen that when the number of particles increases to 500, the volatility of the estimated curves of the PF and IMFA-PF algorithms decreases, and the results of the PF algorithm converge, indicating that the increase in the number of particles can improve the particle degradation phenomenon to a certain extent. In addition, the positioning accuracy and stability of the PF and IMFA-PF algorithms increase with the increase in the number of particles, but IMFA-PF can exceed the positioning effect of the traditional PF algorithm using 500 particles with only 50 particles, effectively reducing the complexity of the algorithm.

[0197] In order to further study the advantages of IMFA-PF algorithm in locating multiple targets, the geographical locations of the five targets observed by the UAV were calculated using PF, BA-PF, SFA-FA and IMFA-PF algorithms. N The average RMSE and average operation time of 10 calculation results are set to 500 for comparative analysis, and the results are shown in Table 3. The target positioning effects of BA-PF, SFA-FA and IMFA-PF are generally better than PF, but the running time increases due to the increase in the computational complexity of the intelligent optimization algorithm. Secondly, as the distance between the target and the imaging center increases, the positioning error also increases. The reason is that the distortion of the pod lens causes the nonlinearity of the observation equation to intensify, and IMFA-PF can improve the particle diversity through multiple mutation operations, while improving the positioning accuracy and alleviating the impact of distortion. Therefore, compared with the BA-PF and SFA-FA algorithms, the positioning error of IMFA-PF at T4 and T5, which are far from the imaging center, is smaller than that of the closer targets, and the positioning accuracy of multiple targets is higher.

[0198] Table 3 Positioning accuracy and operation time of different algorithms

[0199]

[0200] 4 Conclusion

[0201] Aiming at the problem that the UAV target positioning system has serious nonlinear disturbances, resulting in low positioning accuracy, a UAV target positioning method based on improved firefly optimized particle filter is proposed. By combining the attraction and movement mechanism of particle filter and firefly algorithm, the global search ability of particles can be enhanced, so as to effectively adapt to the changes of objective function and the high nonlinearity of positioning system. Furthermore, multiple mutation strategies and elastic force mechanisms are introduced to change the interaction mode of particles, prevent unreasonable aggregation of particles and low-quality convergence caused by over-optimization. Simulation and actual experimental results show that the proposed algorithm can improve the influence of nonlinear disturbance on particle distribution and effectively solve the problem of particle degradation. Compared with existing algorithms, it has better robustness and positioning accuracy.

[0202] The above is an example of the best implementation of the present invention, and the parts not described in detail are common knowledge of ordinary technicians in the field. The protection scope of the present invention shall be based on the content of the claims, and any equivalent transformation based on the technical enlightenment of the present invention is also within the protection scope of the present invention.

Claims

1. UAV target positioning method based on improved mutant firefly optimized particle filter, It is characterized in that The following steps are involved: Step 1: Collect the real-time position observation value of the target, add a random noise to the real-time observation value each time, repeat N times, and obtain N position data. Each position data corresponds to the state of a particle. N particles form a particle swarm. The position state corresponding to the particle swarm at time k-1 constitutes the particle state set Step 2: According to the particle state set at time k-1 Using the state equation of the positioning system to predict the particle state set at time k Calculate the luminous brightness of each particle and select the particle with the largest brightness as the optimal particle. The state value corresponding to the optimal particle is Step 3: Calculate the optimal particle Attraction to other particles The position of each particle at time k is updated and calculated according to the attraction degree. The updated position of the i-th particle is Repeat m for the updated particle swarm 1 The position is swapped and mutated. The particle with the best position is retained each time, and the particles after swapping and mutation are obtained. The state set corresponding to the particles after swapping and mutation is: Step 4: Repeat the m-step for the exchanged and mutated particle swarm. 2 The position differential mutation is performed again to form a new particle swarm after differential mutation. The N particles with the best performance are selected from the particle swarm before the position differential mutation and the new particle swarm after differential mutation as the particle swarm for the next iteration. The state set corresponding to the particle swarm is Step 5: Calculate the density ρ of the particle group obtained in the previous step k , judge ρ k Is it greater than the preset threshold ρ th , then the particle group is dispersed according to the elastic mechanism to obtain the dispersed particle group. The position state of the i-th dispersed particle is Otherwise, keep the position of the particle group unchanged; Step 6: Use the target observation value z at the current k moment k Calculate normalized particle weights Step 7: Output the target location estimate at time k 2. According to claim 1, the unmanned aerial vehicle target positioning method based on improved variant firefly optimized particle filtering, It is characterized in that The formula for updating the position of each particle at time k according to the attraction in step 3 is: In the formula, is the updated state value of particle i after moving at time k; is the initial state value of particle i at time k, α∈[0,1] is the step factor, and rand is a random number uniformly distributed on [0,1]; The calculation formula of α is Where c is the scale factor and μ is the position factor.

3. The unmanned aerial vehicle target positioning method based on improved variant firefly optimized particle filtering according to claim 1, It is characterized in that The specific method of position exchange mutation in step 3 is: A particle is randomly selected and the two dimensions of its state space are exchanged to generate a new particle. The particle with the best position is retained in the new particle and the original particle. The position exchange mutation is completed once, and the mutation operation is repeated m times. 1 Second-rate.

4. The unmanned aerial vehicle target positioning method based on improved variant firefly optimized particle filtering according to claim 1, It is characterized in that The specific method of each position difference mutation in step 4 is: Randomly select two particles to perform differential mutation to obtain new particles. and Perform position differential mutation to obtain new particles The calculation formula is: Where F is the scaling factor.

5. The unmanned aerial vehicle target positioning method based on improved variant firefly optimized particle filtering according to claim 1, It is characterized in that In step 5, ρ k The calculation formula is: Where num is the operation for calculating the number of particles in the judgment area; r j is the radius of the judgment area; N is the total number of particles; The calculation formula is: In the formula, ξ is the rebound coefficient, k is the elastic coefficient, For elasticity, is the compression amount of particle i at time k.

Images