An unmanned aerial vehicle path deception jamming method based on intelligent optimization algorithm

Abstract

The application discloses a method for unmanned aerial vehicle path deception jamming based on an intelligent optimization algorithm, belongs to the technical field of unmanned aerial vehicle path deception jamming, and is used for unmanned aerial vehicle path deception jamming. The method comprises the following steps: establishing a dynamic equation and a kinematic equation of an unmanned aerial vehicle to obtain an expected flight trajectory of the unmanned aerial vehicle; introducing a particle filtering method based on a glowworm swarm optimization algorithm to track and estimate the state of the unmanned aerial vehicle; and introducing an improved particle swarm optimization algorithm to calculate a public time delay in a retransmission type jamming and to interfere with and guide the positioning information of the unmanned aerial vehicle. Compared with the prior art, the application simulates the real flight process and the maneuvering change of the unmanned aerial vehicle, reflects the attitude change and the dynamic response of the unmanned aerial vehicle, alleviates the particle degradation and sample impoverishment problems, improves the tracking stability and the precision under noisy observation, improves the precision of deception positioning and the system robustness, enhances the overall deception jamming effect, and provides technical support for improving the deception jamming capability of an end defense system.

Description

TECHNICAL FIELD

[0001] The application discloses an unmanned aerial vehicle path deception jamming method based on an intelligent optimization algorithm and belongs to the technical field of unmanned aerial vehicle path deception jamming. BACKGROUND

[0002] The core of target tracking technology lies in effective estimation of target states, wherein, as a classical method for optimally estimating states of a linear Gaussian system, Kalman Filter (KF) has been widely applied to information fusion and tracking problems of multi-source sensors such as radars and infrareds due to its minimum mean square error criterion and recursive calculation capability. However, in actual battlefield environments, target movement usually has nonlinear and variable characteristics, and is affected by external interference and sensor measurement errors, so that the traditional KF method has obvious deficiencies in modeling capability and robustness, and is difficult to meet the accurate tracking requirements of modern fire control systems for complex dynamic targets.

[0003] With the wide application of satellite navigation technology, the global satellite navigation positioning system has become an indispensable positioning means in modern civilian and military systems. However, in a complex electromagnetic environment or confrontation scene, satellite signals are easily affected by deception jamming, especially the low-technology threshold and strong concealment of the Meaconing deception jamming, which has attracted widespread attention. This jamming method receives real satellite signals and delays the forwarding, so that the unmanned aerial vehicle cannot distinguish the original and fake signals, thereby misjudging the propagation time, causing the pseudo-range calculation deviation and generating false positioning information. In the actual flight process, due to the noise in the radar observation data, directly using the original measurement trajectory for processing will cause the target state estimation accuracy to be reduced, and affect the implementation of the subsequent path planning and deception jamming strategy. The traditional particle filter method often faces the problems of particle degeneration and sample impoverishment when dealing with nonlinear and non-Gaussian systems, resulting in unstable filtering effect. In addition, in the process of deception jamming, the delay time of the forwarding device is difficult to accurately calculate due to the inability to obtain the real position information of the unmanned aerial vehicle, which affects the accuracy of the deception positioning. In the path jamming strategy, the current algorithm generally has weak global search ability and is easy to fall into local optimum, and it is difficult to generate a deception path with sufficient deception and jamming effect. The Meaconing method does not need to reconstruct the navigation signal, has the advantages of simple implementation and significant jamming effect, and is suitable for the demand of attracting and path deception of target position. SUMMARY

[0004] The application aims to provide an unmanned aerial vehicle path deception jamming method based on an intelligent optimization algorithm to solve the problems of weak global search ability, easy to fall into local optimum and poor jamming effect of the unmanned aerial vehicle path deception jamming algorithm in the prior art.

[0005] An unmanned aerial vehicle path deception jamming method based on an intelligent optimization algorithm comprises:

[0006] S1 establishes the dynamics equation and kinematics equation of the unmanned aerial vehicle, and obtains the expected flight trajectory of the unmanned aerial vehicle;

[0007] S2 introduces an improved particle swarm optimization algorithm to track and estimate the state of the unmanned aerial vehicle;

[0008] S3 introduces a particle filtering method based on glowworm swarm optimization to calculate the public time delay in the relay interference and to interfere and guide the positioning information of the unmanned aerial vehicle.

[0009] S1 includes S1.1, the dynamics of the unmanned aerial vehicle is analyzed as a rigid body, and the mass center dynamics equation of the unmanned aerial vehicle in the body coordinate system is:

[0010] ;

[0011] In the formula, 、 and are the roll angular velocity, yaw angular velocity and pitch angular velocity respectively, 、 and are the velocity components of the velocity vector in the axis, axis and axis respectively, is the engine thrust, 、 and are the three direction force drag, lift and side force decomposed by the air force in the airflow coordinate system, is the angle between the thrust vector and the longitudinal axis of the body, is the mass of the unmanned aerial vehicle, is the gravitational acceleration, is time, is the pitch angle, is the roll angle, is the angle of attack, is the side slip angle;

[0012] S1.2, the dynamics equation of the unmanned aerial vehicle rotating around the mass center in the body coordinate system is:

[0013] ;

[0014] ;

[0015] In the formula, 、 、 are the moments of inertia, is the inertia product, is the variable that can be solved, , and are the external moment acting on the UAV in the body coordinate system axis, axis and axis.

[0016] S1 includes S1.3 to establish the relationship between the velocity components of the UAV in the ground coordinate system:

[0017] ;

[0018] wherein, , and are the three-dimensional coordinates of the UAV in axis, axis and axis, , , , , and are the three-dimensional velocities of the UAV in axis, axis and axis, is the yaw angle;

[0019] S1.4 establishes the relationship between the angular velocity components of the UAV in the body coordinate system, which is equivalent to the kinematic equation of the UAV rotating around the center of mass:

[0020] ;

[0021] wherein, , are the derivatives of the roll angular velocity, yaw angular velocity and pitch angular velocity, respectively, and the state quantity of the UAV is , is the control input.

[0022] S2 includes S2.1 to add observation errors to the position information matrix at time as the observation value of the state transition process :

[0023] = ;

[0024] wherein, , , , , and represent the added observation error coefficients, a random number representing a normal distribution with mean 0 and variance 1;

[0025] Calculate the particle filter;

[0026] At the initial time, sample the system state to generate an initial particle set, pass the prior distribution of the system Initialize the particles and assign equal initial weights to each particle:

[0027] ;

[0028] where, is the state of the i-th particle at the initial time, is the total number of particles, is the initial state, is the weight of the i-th particle at the initial time; At the time t, sample each particle according to the state transition model, and pass the particle distribution at the time t to the current time:

[0029] ;

[0030] ;

[0031] where, is the state of the i-th particle at the time t, is the particle distribution at the time t passed to the prior distribution at the time t, is the state value at the time t; Obtain the conditional posterior probability density of the state characterized by the weighted particles : ;

[0032]

[0033] ;

[0034] where, is the observation value at the time t, is the weight of the i-th particle at the time t, is the Dirac function; Adjust the weights of the particles using the current observation data:

[0035] ;

[0036] ;

[0037] where,​​​​​​​​​​ is the proportionality symbol, is the particle state is the observation data is the probability;

[0038] The updated particle weights are normalized:

[0039] ;

[0040] where, is the normalized weight, is the i-th particle;

[0041] The state estimate of the target is the weighted average of the current state distribution:

[0042] ;

[0043] where, is the state estimate at time t;

[0044] Finally, the error covariance is calculated:

[0045] .

[0046] S2 includes S2.2 calculating the attraction degree value between each particle and the global optimal firefly individual, updating the position of the particle;

[0047] In the brightness calculation of the firefly algorithm, the of the particle filter is introduced, and the higher the brightness value, the better the spatial position, and the greater the weight of the corresponding particle, and it is closer to the true position;

[0048] The brightness of the firefly individual is calculated as:

[0049] ;

[0050] where, is the observation noise covariance matrix;

[0051] According to the brightness value, the calculation of the attraction degree of the firefly individual is:

[0052] ;

[0053] where, is the maximum attraction degree value, taking the value of 1; is the light attenuation coefficient, taking the value range of ; represents the distance between two fireflies,​​​ It is a natural constant;

[0054] exist At any moment, assuming For individuals with lower brightness values, For individuals with higher brightness values, through The globally optimal firefly individual at time t is redefined as follows:

[0055] ;

[0056] In the formula, for time Location, This is the step size factor, and its value range is... , It is a random number.

[0057] S2 includes S2.3. After updating the particle positions, the brightness value of the particle is calculated and compared with the current global optimum, dynamically updating the optimal solution. If the global optimum meets the set convergence threshold, the optimal particle is selected as the final state estimation result; otherwise, the particle swarm is resampled. The resampling process is as follows:

[0058] Construct the cumulative distribution function based on the normalized weights. ,make , ;

[0059] ;

[0060] Particles are selected through random sampling. Generate within the interval random numbers The cumulative distribution function is used to map random numbers. For each random number... Find the index of the cumulative distribution function. ,satisfy:

[0061] ;

[0062] Select the individual particles As a particle in the new particle set, it continues to perform the update process based on the firefly algorithm in the new particle set until the convergence condition is met.

[0063] S3 includes the S3.1 drone and the... The actual spatial distance between the satellites is The drone measured the same as the first The distance between the satellites is , The relationship between the The relationship between the

[0064]

[0065] In the formula, The coordinates of the first The position information in the final output of the improved particle swarm optimization algorithm is taken as the coordinates of the unmanned aerial vehicle , , The propagation delay error caused by the ionosphere and troposphere of the atmosphere, The clock error of the unmanned aerial vehicle, The clock error of the satellite, The speed of light;

[0066] In the repeater deception jamming, the signal delay is artificially increased, the pseudo-range measurement value of the unmanned aerial vehicle is changed, and the relationship between the The relationship between the The relationship between the

[0067]

[0068] In the formula, The signal delay amount of the first satellite introduced;

[0069] The repeater coordinates After the satellite signal is received by the repeater device, the satellite position is calculated, and the propagation time delay of the repeater signal is obtained according to the satellite position and the repeater coordinates According to the satellite position coordinates and the position coordinates of the repeater device, the propagation time delay of the satellite signal to the repeater device is calculated According to the actual position of the unmanned aerial vehicle and the position of the repeater device, the propagation time delay of the signal is calculated The actual position of the unmanned aerial vehicle is obtained by the detection device;

[0070] The repeater time delay amount of each satellite is:

[0071]

[0072] According to the geometric relationship of the satellite, the repeater device and the unmanned aerial vehicle position, the common time delay amount is introduced, which satisfies , The total number of the repeater satellite signals is:

[0073] After the common time delay amount is introduced, the repeater pseudo-range measurement equation is:

[0074] ​​​​

[0075] According to the principle of satellite positioning solution, the introduction of common time delay causes the change of positioning solution clock error, and the clock error of the unmanned aerial vehicle becomes:

[0076] ;

[0077] In the formula, is the common time delay after the introduction .

[0078] S3 includes S3.2 optimizing the common time delay with the objective function :

[0079] ;

[0080] In the formula, is the objective function, is the distance between the preset deception positioning position of the jammer and the actual deception positioning position, is the preset deception positioning position of the jammer, is the actual deception positioning position of the jammer, is the decision vector composed of the common time delay , and satisfies , is the corrected forwarding time delay.

[0081] S3 includes S3.3 seeking the optimal solution by using an intelligent optimization algorithm based on particle filtering and Levy flight strategy;

[0082] S3.3.1 initializes the algorithm parameters, sets the maximum size of the population of the algorithm , the maximum number of iterations , the upper bound and lower bound of the search range , , the number of hyperparameters of the optimization problem ; the mutation probability in genetic algorithm ; the inertia weight , learning factor and , random factor and are random numbers in the interval , particle velocity ; the control parameter in Levy flight strategy ; the mapping coefficient in chaos mapping strategy ;

[0083] S3.3.2 initializes all population individual parameters based on the Logistic chaos mapping function and the reference individual in the real number coding mode ​, the expression of Logistic function is:

[0084] ;

[0085] where, is the value of the th individual;

[0086] The individual beyond the upper and lower bounds of the search range is defined as the average value of the upper and lower bounds;

[0087] S3.3.3 Calculate the fitness function value of all individuals according to the function, the fitness function constructs the distance between the preset deception positioning position of the jammer and the actual deception positioning position, save the individual optimal position in the iteration process and the global optimal position ;

[0088] S3.3.4 According to the set , select part of the population individuals, introduce the Levy flight strategy for mutation operation:

[0089] ;

[0090] where, represents the current iteration number, is the th particle in the th iteration position, is the scaling factor, which controls the Levy flight jump amplitude, , is a random vector of Levy flight distribution, is the tensor product operator symbol, and are control parameters, which change with the iteration number, the expression is:

[0091] ;

[0092] ;

[0093] where, is the maximum iteration number of the algorithm;

[0094] For the remaining part of the population individuals, continue to update the speed and position through the particle swarm algorithm:

[0095] ;

[0096] where, is the th particle in the th iteration speed; controls the influence of historical speed; and Adjust the weight of individual and group experience; and Enhance search randomness;

[0097] The position updating formula of the particle updates the current position of the particle by moving the particle in the current speed direction by a set step size:

[0098] ;

[0099] S3.3.5 Calculate the updated fitness value of all individuals, the individual optimal position and the global optimal position;

[0100] S3.3.6 If the maximum number of iterations is reached or the threshold condition is met, the algorithm ends and the common time delay is obtained ; Otherwise, return to S3.3.4 for the next iteration loop.

[0101] The common time delay is obtained as input to obtain the clock error of the unmanned aerial vehicle and the that lures the unmanned aerial vehicle from point A to point B, so that the unmanned aerial vehicle obtains false pseudorange information, and the positioning position calculated by the unmanned aerial vehicle is point B, thereby realizing the relay deception jamming of the unmanned aerial vehicle.

[0102] Compared with the prior art, the present application has the following beneficial effects: the six-degree-of-freedom model of the unmanned aerial vehicle is modeled and simulated, so that the real flight process and maneuvering changes of the unmanned aerial vehicle are more accurately simulated, and the model can better reflect the attitude changes and dynamic responses of the unmanned aerial vehicle in flight; the particle filtering method based on glowworm optimization is introduced, the particle sampling strategy and weight updating mechanism are improved, the problems of particle degradation and sample impoverishment are effectively alleviated, and the tracking stability and precision under noisy observation are significantly improved; the common time delay optimization calculation model based on the improved particle swarm optimization algorithm is constructed, the relay error is dynamically compensated and corrected, thereby improving the precision of deception positioning and the system robustness, and enhancing the overall deception jamming effect; the system is constructed to integrate trajectory generation, state tracking and path deception into an integrated unmanned aerial vehicle deception guidance, and provides effective technical support for improving the deception jamming capability of the terminal defense system. BRIEF DESCRIPTION OF DRAWINGS

[0103] Figure 1 is a three-dimensional flight trajectory diagram of the unmanned aerial vehicle;

[0104] Figure 2 is a curve diagram of the axis position change of the unmanned aerial vehicle in the flight process ;

[0105] Figure 3 is a curve diagram of the axis position change of the unmanned aerial vehicle in the flight process The axis position change curve diagram;

[0106] Figure 4 is the change curve diagram of the angle of attack and the sideslip angle during the flight of the unmanned aerial vehicle The axis position change curve diagram;

[0107] Figure 5 is the change curve diagram of the angle of attack and the sideslip angle during the flight of the unmanned aerial vehicle

[0108] Figure 6 is the change curve diagram of the pitch angle, the yaw angle and the roll angle during the flight of the unmanned aerial vehicle

[0109] Figure 7 is the change curve diagram of the roll angle velocity, the yaw angle velocity and the pitch angle velocity during the flight of the unmanned aerial vehicle

[0110] Figure 8 is the change curve diagram of the inertial position during the flight of the unmanned aerial vehicle

[0111] Figure 9 is the Y-axis error diagram before filtering

[0112] Figure 10 is the Y-axis error diagram after filtering

[0113] Figure 11 is the comparison diagram of the three-axis errors before and after filtering

[0114] Figure 12 is the three-dimensional flight trajectory diagram of the unmanned aerial vehicle under the retransmission type deception jamming. DETAILED DESCRIPTION

[0115] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the present application are described clearly and completely below. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the present application.

[0116] A method for deception jamming of an unmanned aerial vehicle path based on an intelligent optimization algorithm, comprising:

[0117] S1, establishing a dynamic equation and a kinematic equation of the unmanned aerial vehicle, and obtaining a desired flight trajectory of the unmanned aerial vehicle;

[0118] S2, introducing an improved particle swarm optimization algorithm to track and estimate the state of the unmanned aerial vehicle;

[0119] S3, introducing a particle filtering method based on the glowworm swarm optimization to calculate the public time delay in the retransmission type jamming, and interfering and guiding the positioning information of the unmanned aerial vehicle.

[0120] S1 includes S1.1 analyzing the dynamics of the unmanned aerial vehicle as a rigid body, the mass center dynamics equation of the unmanned aerial vehicle in the body coordinate system is:

[0121] ;

[0122] wherein, , and are roll angular velocity, yaw angular velocity and pitch angular velocity respectively, , and are the velocity components of the velocity vector in the axis, axis and axis, is the engine thrust, , and are the force drag, lift and side force in three directions decomposed in the airflow coordinate system, is the angle between the thrust vector and the body longitudinal axis, is the mass of the unmanned aerial vehicle, is the acceleration of gravity, is time, is the pitch angle, is the roll angle, is the angle of attack, is the angle of sideslip;

[0123] S1.2 the dynamics equation of the unmanned aerial vehicle rotating around the mass center in the body coordinate system is:

[0124] ;

[0125] ;

[0126] wherein, , , are the moments of inertia, is the inertia product, is the variable to be solved, , and are the components of the external moment acting on the unmanned aerial vehicle in the body coordinate system in the axis, axis and axis.

[0127] S1 includes S1.3 establishing the relationship between the velocity components of the unmanned aerial vehicle in the ground coordinate system:

[0128] ;

[0129] wherein, , and are three-dimensional coordinates of the UAV in axis, axis and axis, , , , , and are three-dimensional velocities of the UAV in axis, axis and axis, is the yaw angle;

[0130] S1.4 establishes the relationship between the angular velocity components of the UAV in the body coordinate system, which is equivalent to the kinematic equation of the UAV rotating around the center of mass:

[0131] ;

[0132] wherein, , are the derivatives of the roll angle velocity, the yaw angle velocity and the pitch angle velocity, respectively, and the state quantity of the UAV is , is the control input.

[0133] S2 includes S2.1 adding observation errors to the position information matrix at the time as the observation value of the state transition process :

[0134] = ;

[0135] wherein, , , , , and are the added observation error coefficients, is a random number of a normal distribution with a mean of 0 and a variance of 1;

[0136] calculating the particle filter;

[0137] At the initial time, the system state is sampled to generate an initial particle set, the particles are initialized through the prior distribution of the system , and each particle is assigned an equal initial weight:

[0138] ;

[0139] In the formula, For the initial time, the first The state of each particle The total number of particles, This is the initial state. For the initial time, the first The weight of each particle;

[0140] exist At each moment, according to the state transition model, each particle is sampled, and The particle distribution at time 1 is passed to the current time:

[0141] ;

[0142] In the formula, for Time of the first The state of each particle for The particle distribution at time is transmitted to Prior distribution at time, for The state value at any given time;

[0143] Obtain the conditional posterior probability density of the state using weighted particles. :

[0144] ;

[0145] In the formula, for The observed value at time, for Time of the first The weight of each particle, It is the Dirac function;

[0146] Using current observation data Adjust the particle weights:

[0147] ;

[0148] In the formula, It is a proportional symbol. In particle state Time observation data The probability of;

[0149] Normalize the updated particle weights:

[0150] ;

[0151] In the formula, The normalized weights For the first One particle;

[0152] The target's state estimate is the weighted average of the current state distribution:

[0153] ;

[0154] In the formula, for State estimate at time 1;

[0155] Finally, the error covariance was calculated. :

[0156] .

[0157] S2 includes S2.2, which calculates the attraction value between each particle and the globally optimal firefly individual, and updates the particle's position;

[0158] Particle filtering is introduced into the brightness calculation of the firefly algorithm. A higher brightness value indicates a better spatial position, a greater weight for the corresponding particle, and a closer resemblance to its actual position.

[0159] The brightness of an individual firefly The calculation is as follows:

[0160] ;

[0161] In the formula, To observe the noise covariance matrix;

[0162] Based on brightness values, the attractiveness of individual fireflies The calculation is as follows:

[0163] ;

[0164] In the formula, The maximum attraction value is set to 1. The optical attenuation coefficient has a value range of [value range missing]. ; This indicates the distance between two fireflies. It is a natural constant;

[0165] exist At any moment, assuming For individuals with lower brightness values, For individuals with higher brightness values, through The globally optimal firefly individual at time t is redefined as follows:

[0166] ;

[0167] In the formula, For the moment of the position, is a step factor, the value range is , is a random number.

[0168] S2 includes S2.3 updates the particle position, calculates the brightness value of the particle, and compares with the current global optimal value, dynamically updates the optimal solution; if the global optimal value meets the set convergence threshold, the optimal particle is selected as the final state estimation result, otherwise, the particle swarm is resampled, and the resampling process is:

[0169] According to the normalized weight, the cumulative distribution function is constructed , , ;

[0170] ;

[0171] Select particles by random sampling, generate random numbers in the interval , map the random numbers using the cumulative distribution function, for each random number , find the index of the cumulative distribution function, which satisfies:

[0172] ;

[0173] Select the th particle as a particle in the new particle set, continue to execute the update process based on the glowworm algorithm in the new particle set until the convergence condition is met.

[0174] S3 includes S3.1 the real space distance between the unmanned aerial vehicle and the first satellite is , the distance measured by the unmanned aerial vehicle and the first satellite is , The relationship between , and is:

[0175] ;

[0176] In the formula, is the coordinate of the first satellite, and the position information in output by the improved particle swarm optimization algorithm is taken as the coordinate of the unmanned aerial vehicle , The propagation delay error caused by the atmospheric ionosphere troposphere, The clock error of the UAV, The clock error of the satellite, The speed of light;

[0177] In the repeater spoofing jamming, the signal delay is artificially increased, the pseudo-range measurement value of the UAV is changed, and the relationship between and is calculated as follows:

[0178] ;

[0179] In the formula, is the signal delay amount of the first satellite introduced;

[0180] The repeater coordinates are determined After the satellite signal is received by the repeater, the satellite position is calculated, and the propagation time delay of the repeater signal is obtained according to the satellite position and the repeater coordinates ; according to the satellite position coordinates and the position coordinates of the repeater, the propagation time delay of the satellite signal reaching the repeater is calculated ; according to the actual position of the UAV and the position of the repeater, the propagation time delay of the signal is calculated

[0181] The repeater time delay amount of each satellite is:

[0182] ;

[0183] According to the geometric relationship of the satellite, the repeater, and the UAV positions, a common time delay amount is introduced, which satisfies , is the total number of satellite signals that are repeated;

[0184] After the common time delay amount is introduced, the repeater pseudo-range measurement equation is:

[0185] ;

[0186] According to the satellite positioning solution principle, the introduction of the common time delay amount causes the change of the positioning solution clock error, and the clock error of the UAV becomes:

[0187] ;

[0188] In the formula, is the after the common time delay amount is introduced.

[0189] S3 includes S3.2 using the objective function Optimize common latency :

[0190] ;

[0191] In the formula, Let be the objective function, and be the distance between the preset deception positioning position and the actual deception positioning position of the jammer. The preset deception positioning location for the jammer. The actual deceptive location of the jammer. It is determined by the common delay. The constructed decision vector satisfies , This is the corrected forwarding delay.

[0192] S3, including S3.3, employs an intelligent optimization algorithm based on particle filtering and the Levy flight strategy to seek the optimal solution;

[0193] S3.3.1 Initialize algorithm parameters and set the maximum population size for the algorithm. Maximum number of iterations Upper bound of the search range and the lower world Number of hyperparameters in the optimization problem Mutation probability in genetic algorithms Inertia weights in particle swarm optimization Learning factors and Random factors and yes Random numbers in the interval, particle velocity Control parameters in Levi's flight strategy Mapping coefficients in chaotic mapping strategies ;

[0194] S3.3.2 uses real-number encoding to initialize all population individual parameters based on the Logistic chaotic mapping function and the baseline individual. The expression for the Logistic function is:

[0195] ;

[0196] In the formula, For the first Individual values;

[0197] Individuals that fall outside the upper and lower bounds of the search range are defined as the average of the upper and lower bounds;

[0198] S3.3.3 According to The function calculates the fitness function value for all individuals. The fitness function is used to construct the distance between the deception location preset by the interference machine and the actual deception location. The optimal position of each individual during the iteration process is saved. and the global optimal position ;

[0199] S3.3.4 According to the settings Select a subset of individuals from the population and introduce the Levy flight strategy for mutation manipulation:

[0200] ;

[0201] In the formula, Indicates the current iteration number. For the first The particle in the first The position of the next iteration. The scaling factor controls the jump amplitude of Levi's flight. , Let Lévy's flight distribution be a random vector. This is the tensor product operator. and To control the parameters, which change with the number of iterations, the expression is:

[0202] ;

[0203] ;

[0204] In the formula, This represents the maximum number of iterations for the algorithm.

[0205] For the remaining individuals in the population, continue to update their velocity and position using the particle swarm optimization algorithm:

[0206] ;

[0207] In the formula, It is the first The particle in the first The speed of each iteration; Controlling the impact of historical speed; and Adjusting the weight of individual and group experiences; and Enhance search randomness;

[0208] The particle position update formula updates the particle's current position by moving the particle a set step size in the direction of its current velocity:

[0209] ;

[0210] S3.3.5 Calculate all individual updated fitness values, individual optimal positions and global optimal position;

[0211] S3.3.6 If the maximum iteration number is reached or the threshold condition is met, the algorithm ends and the common time delay is obtained ; Otherwise, return to S3.3.4 for the next iteration loop.

[0212] The common time delay is obtained , as input to obtain the clock error of the unmanned aerial vehicle and the , so that the unmanned aerial vehicle obtains false pseudorange information, and the positioning position calculated by the unmanned aerial vehicle is B point, so as to realize the relay deception jamming of the unmanned aerial vehicle.

[0213] In actual flight process, the radar system inevitably introduces noise and error when measuring the trajectory of the unmanned aerial vehicle. Therefore, it is necessary to superimpose an observation error model on the ideal trajectory to construct observation data with certain error. For this observation data, the improved particle swarm optimization algorithm and the particle filtering method based on firefly algorithm optimization are introduced to cooperatively improve the accuracy and robustness of the state estimation of the unmanned aerial vehicle, so as to realize accurate tracking and error suppression of the target trajectory. The high-precision estimated trajectory output by the above particle filtering system is used as the deception jamming input, and the common time delay in the interference signal is optimized and calculated by the particle swarm optimization algorithm. With the help of the method, the target can be misled to the preset error position, so that after the position data with error is received by the opponent system, significant positioning deviation is caused, and the purpose of effective interference and tactical deception is achieved.

[0214] Based on the established unmanned aerial vehicle dynamics and kinematics model, any flight trajectory can be generated to verify the adaptability and effectiveness of the model. As shown in Figure 1 , the following flight path is preset: the unmanned aerial vehicle first flies horizontally along the axis for 50 seconds, then climbs to a height of 50 seconds, and finally hovers within a certain range. The position change curves of the axis, axis and axis during the flight of the unmanned aerial vehicle are shown in Figure 2 , Figure 3 and Figure 4 . The flight process simulates the complete task flow of the unmanned aerial vehicle from take-off, flight to the target area and implementation of reconnaissance.

[0215] Based on the established model, real-time monitoring and recording of each state variable of the unmanned aerial vehicle during flight can be realized. The curves of the angle of attack and the side slip angle changing with time are shown in Figure 5 , alpha represents the angle of attack, and beta represents the side slip angle; the curves of the pitch angle, yaw angle and roll angle changing with time are shown in Figure 6As shown, phi represents the pitch angle, theta represents the yaw angle, and psi represents the roll angle; the roll angle velocity, the yaw angle velocity, and the pitch angle velocity curves over time are as shown in Figure 7 As shown, p represents the roll angle velocity, q represents the yaw angle velocity, and r represents the pitch angle velocity; the inertial position curve over time is as shown in Figure 8 As shown, Xe represents the axis position, Ye represents the axis position, and Ze represents the axis position; the simulation results show that the change trend of each state quantity conforms to the basic laws of aerodynamics and flight dynamics, has good physical consistency, can truly reflect the dynamic characteristics of the unmanned aerial vehicle in actual flight, and verifies the rationality and effectiveness of the model.

[0216] To verify the effectiveness of the proposed particle filtering method based on glowworm algorithm optimization in the state estimation of the unmanned aerial vehicle, a comparative analysis is performed in terms of trajectory error estimation. Figure 9 , Figure 10 The error in the X-axis direction before and after the particle filtering method based on glowworm algorithm optimization is respectively shown in the error over time. Before filtering, the error gradually increases over time, and the fluctuation is significant, indicating that there is a large observation noise and estimation bias in the system without filtering. After using the PF method, the trajectory error is greatly reduced, and the error interval is obviously contracted, indicating that the PF has certain filtering ability and robustness.

[0217] Further, the comparison chart of the error in the three-axis directions before and after filtering is given in Figure 11 . The blue points in the chart represent the error in the unfiltered state, and the red points represent the error after the PF method. As can be seen from the chart, the PF method effectively reduces the fluctuation range of the error, making the error distribution more concentrated. This is because the brightness function based on the observation value and the attraction degree mechanism are introduced in the particle updating process, so that the particles can adaptively tend to the estimation direction with smaller error, significantly improving the estimation accuracy.

[0218] Based on the relay satellite signal deception model, the flight trajectory of the unmanned aerial vehicle in the interference environment is simulated and visualized using MATLAB. The initial position, target point, and coordinates of four satellites are set, the signal propagation delay is optimized by the improved particle swarm optimization algorithm, the navigation system is misled, the unmanned aerial vehicle gradually deviates from the real trajectory in control, and the three-dimensional flight trajectory of the unmanned aerial vehicle under relay deception interference is as shown in Figure 12 .

[0219] The above examples are only used for illustrating the technical solutions of the present application, and are not intended to limit the present application. Although the present application has been described in detail with reference to the foregoing examples, it should be understood by those skilled in the art that the technical solutions recorded in the foregoing examples can be modified, or some or all of the technical features can be replaced by equivalent replacements, and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for unmanned aerial vehicle path deception jamming based on intelligent optimization algorithm, characterized in that, Comprise: S1 establishes the dynamic equation and kinematics equation of unmanned aerial vehicle, obtains the expected flight trajectory of unmanned aerial vehicle; S2 introduces improved particle swarm optimization algorithm, tracks and estimates the state of unmanned aerial vehicle; S3 introduces particle filtering method based on glowworm algorithm optimization to calculate the public time delay in the relay interference, and interferes and guides the positioning information of unmanned aerial vehicle; S2 includes S2.1 calculates error covariance; S2.2 calculates the attraction value between each particle and the global optimal glowworm individual, and updates the position of the particle; S2.3 after updating the position of the particle, the brightness value of the particle is calculated and compared with the current global optimal value, and the optimal solution is dynamically updated; if the global optimal value meets the set convergence threshold, the optimal particle is selected as the final state estimation result, otherwise, the particle swarm is resampled; S3 includes S3.1 calculates the clock difference of unmanned aerial vehicle after introducing public time delay; S3.2 Objective function Optimizing common latency ; S3.3 adopts intelligent optimization algorithm based on particle filtering and Levy flight strategy to seek optimal solution; S3.3.1 initialize algorithm parameters, set the maximum size of the population of the algorithm , the maximum number of iterations , the upper bound of the search range and the lower bound , the number of hyperparameters of the optimization problem ; mutation probability in genetic algorithm ; inertia weight in particle swarm optimization algorithm , learning factor and , random factor and are random numbers in the interval, particle velocity ; Control parameters in Levy flight strategy Mapping coefficients in chaotic mapping strategy ; S3.3.2 All population individual parameters are initialized based on the Logistic chaotic mapping function with the real number coding The expression of the Logistic function is: ； In the formula, is the first individual value; The individual exceeding the upper and lower bounds of the search range is defined as the average value of the upper and lower bounds; S3.3.3 according to The function calculates the fitness function value of all individuals, the fitness function constructs the distance between the deception positioning position preset by the jammer and the actual deception positioning position, and saves the individual optimal position in the iteration process and the global optimal position ; S3.3.4 According to the setting Select part of the population individuals, introduce Levy flight strategy for mutation operation: ； wherein, denotes the current iteration number, is the position of the th particle at the th iteration, is a scaling factor that controls the amplitude of the Levy flight jumps, , is a random vector of the Levy flight distribution, is the tensor product operator, and are control parameters that vary with the iteration number, and are expressed as: ； ； In the formula, is the maximum number of iterations for the algorithm; For the remaining individuals in the population, continue to update the speed and position through particle swarm algorithm: ； wherein is the velocity of the th particle at the th iteration; controls the influence of the history velocity; and adjusts the weight of individual and group experience; and enhances search randomness; The position update formula of the particle is to move the current position of the particle by a certain step length in the direction of the current speed: ； S3.3.5 calculates the updated fitness value of all individuals, individual optimal position and global optimal position; S3.3.6 If the maximum number of iterations is reached or the threshold condition is met, the algorithm ends and the common delay quantity is obtained ; otherwise go back to S3.3.4 for the next iteration loop; Obtaining common time delay amount , as input, the clock difference of the UAV and the , making the UAV get the wrong pseudo-range information, and further making the positioning position calculated by the UAV be the B point, so as to realize the retransmission type deception jamming of the UAV.

2. The method of claim 1, wherein the method is based on a smart optimization algorithm. S1 includes S1.1 analyzes the dynamics of unmanned aerial vehicle as a rigid body, and the mass center dynamics equation of unmanned aerial vehicle in body coordinate system is: ； wherein, , and are roll, yaw and pitch angular velocities, respectively, , and are velocity components of the velocity vector in axis, axis and axis, respectively, is the engine thrust, , and are the force drag, lift and side force decomposed in three directions under the airflow coordinate system, is the angle between the thrust vector and the longitudinal axis of the aircraft, is the mass of the UAV, is the gravitational acceleration, is time, is the pitch angle, is the roll angle, is the angle of attack, is the angle of sideslip; S1.2 establishes the dynamics equation of unmanned aerial vehicle rotating around mass center in body coordinate system: ； ； wherein, , , is the moment of inertia, is the product of inertia, is the variable to be solved, , and are the components of the external moment acting on the UAV in the body coordinate system of the axis, axis and axis.

3. The method of claim 2, wherein, S1 includes S1.3 establishes the relationship between velocity components of unmanned aerial vehicle in ground coordinate system: ； wherein, , and are the three-dimensional coordinates of the UAV in the axis, axis and axis respectively, , , , , and are the three-dimensional velocities of the UAV in the axis, axis and axis respectively, is the yaw angle; S1.4 establishes the relationship between angular velocity components of unmanned aerial vehicle in body coordinate system, which is equivalent to the kinematics equation of unmanned aerial vehicle rotating around mass center: ； In the formula, , are the derivatives of roll, yaw and pitch angular velocities, respectively, and the state variables of the UAV are obtained as , is the control input.

4. The method of claim 3, wherein, S2 comprises S2.1 pair matrix of position information at time instant adding observation errors as observations of the state transition process : = ； wherein , , , , and denote the added observation error coefficient, denotes a random number of a normal distribution with mean 0 and variance 1. Calculate particle filtering; At the initial time, the system state is sampled to generate an initial set of particles, by the prior distribution of the system Initialize the particles and assign each particle an equal initial weight: ； wherein is the state of the particle at the initial time step, is the state of the particle at the initial time step, is the total number of particles, is the initial state, is the weight of the particle at the initial time step, is the weight of the particle at the initial time step, At each time step, according to the state transition model, each particle is sampled to propagate the particle distribution from time step t to time step t + 1. At each time step, according to the state transition model, each particle is sampled to propagate the particle distribution from time step t to time step t + 1. At each time step, according to the state transition model, ； wherein is the state of the particle at time is the particle distribution at time is transferred to the prior distribution at time is the state value at time Obtaining a conditional posterior probability density of a weighted particle representation of state : ； wherein is the observation value at time is the weight of the th particle at time is the Dirac function; Utilizing current observation data Adjusting weights of particles: ； wherein is a proportional symbol, is a particle state is observation data is a probability; The updated particle weight is normalized: ； In the formula, is the normalized weight, is the normalized weight, is the normalized weight, The state estimation value of the target is the weighted average value of the current state distribution: ； In the formula, is state estimate at the time instant Finally, the error covariance is calculated : 。 5. The method of claim 4, wherein, S2 comprises S2.2 introducing particle filter in the brightness calculation of glowworm algorithm , the higher the brightness value, the better the spatial position, the greater the weight of the corresponding particle, and the closer to the true position; Brightness of individual fireflies The calculation is: ； In the formula, R is the observation noise covariance matrix; The attractiveness of the firefly individual according to the luminance value is calculated as: ； wherein, is the maximum attraction value, which is equal to 1 ; is the light attenuation coefficient, which is in the range ; denotes the distance between two fireflies, is a natural constant; At the moment, assume for the individual with lower luminance value, for the individual with higher luminance value, redefine the movement process of the individual by the global optimal firefly individual at the moment ； In the formula, is the time of the position, is a step factor, and the value range is , is a random number.

6. The method of claim 5, wherein, S2 includes S2.3 resampling process is: Constructing a cumulative distribution function from normalized weights , such that , ; ； The particles are selected by means of random sampling, in the interval [0, 1] is generated a random number , the random number is mapped using the cumulative distribution function, for each random number , the index of the cumulative distribution function is found , which satisfies: ； selecting a first particle of the first particle As one particle in the new particle set, the new particle set continues to perform the update process based on the glowworm algorithm until the convergence condition is met.

7. The method of claim 6, wherein the method is based on a smart optimization algorithm. S3 includes the S3.1 drone and the... The actual spatial distance between the satellites is The drone measured the same as the first The distance between the satellites is , and The relationship is as follows: ； In the formula, is the first coordinates of the geostationary satellite, the position information in the final output of the improved particle swarm optimization algorithm is taken as the coordinates of the unmanned aerial vehicle , is the propagation delay error caused by the ionosphere and troposphere, is, is the satellite clock error, is the speed of light;​ In the repeater deception jamming, the signal delay is artificially increased, the pseudo-range measurement value of the UAV is changed, and the relationship between the pseudo-range measurement value and the real distance is calculated The relationship between the pseudo-range measurement value and the real distance is calculated The relationship between the pseudo-range measurement value and the real distance is calculated ； In the formula, is the first Signal delay amount of the geostationary satellite; Determine the forwarding coordinates After the satellite signal is received by the forwarding device, the satellite position is calculated, and the propagation delay of the forwarding signal is calculated according to the satellite position and the forwarding coordinates According to the satellite position coordinates and the position coordinates of the forwarding device, the propagation delay of the satellite signal reaching the forwarding device is calculated According to the actual position of the unmanned aerial vehicle and the position of the forwarding device, the propagation delay of the signal is calculated The actual position of the unmanned aerial vehicle is obtained by the detection device The relay time delay of each satellite is: ； According to the geometric relationship of the satellite, the repeater device, and the unmanned aerial vehicle position, a common time delay quantity is introduced , satisfies , is the total number of the satellite signals that are repeated When introducing common time delay After that, the pseudo-range measurement equation is forwarded as: ； According to the principle of satellite positioning solution, the introduction of common time delay causes the change of positioning solution clock difference, and the clock difference of the unmanned plane becomes: ； In the formula, is the public delay amount introduced .

8. The method of claim 7, wherein, S3 includes S3.2: ； In the formula, is the distance between the preset deception positioning position of the jammer and the actual deception positioning position, is the preset deception positioning position of the jammer, is the actual deception positioning position of the jammer, is a decision vector composed of public time delay quantities satisfies , is the modified forwarding time delay.

Images