Unmanned aerial vehicle cluster cooperative target encirclement method based on brain-like computing

Abstract

The application discloses a kind of unmanned aerial vehicle cluster cooperative target encirclement methods based on brain-like computing, belong to autonomous navigation robot and multi-robot coordination control technical field, method includes: obtaining multi-view image under camera visual angle, obtain the local pose and covariance of each unmanned aerial vehicle;According to the local pose and covariance of each unmanned aerial vehicle, obtain the global consistent relative pose of unmanned aerial vehicle cluster;Using unscented Kalman filtering algorithm, the state of target is estimated, and the state information of target in unified coordinate system is obtained;Using Bezier curve generation algorithm, the target motion trajectory under a period of history state is obtained;The method of centroid extension is used to expand the unmanned aerial vehicle cluster, and the encirclement queue of target encirclement is obtained, the least encirclement cost is solved using Gauss Newton method, and the encirclement position of each unmanned aerial vehicle is generated;The motion trajectory of unmanned aerial vehicle is generated using hybrid A* search, and the final encirclement trajectory sequence is generated, which significantly improves the accuracy of motion state estimation.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of autonomous navigation robots and multi-robot coordination control, and particularly relates to a UAV cluster cooperative target encircling method based on brain-like computing. BACKGROUND

[0002] With the deep evolution of information technology, the battlefield situation presents the characteristics of high dynamicity and intense confrontation. Under this background, rapid identification, continuous tracking and accurate capture of high-speed maneuvering targets have become the core ability requirements for improving the efficiency of tactical reconnaissance and building an active defense system. The UAV swarm, with its distributed cooperation advantage, strong spatial coverage capability and flexible formation reorganization characteristics, provides a revolutionary solution for target acquisition tasks in complex battlefield environments. Traditional UAV swarm cooperative capture models are usually based on two key assumptions: first, relying on global absolute position information provided by the global positioning system or similar infrastructure; second, the target's motion pattern is relatively stable or predictable, and the capture process mainly adopts a pre-set fixed geometric formation (such as spherical enclosure, conical interception, etc.). These methods perform well in open areas and low-speed target scenarios. However, when the combat environment evolves into complex urban canyons, jungles or underground spaces, where there is a lack of GPS signal and dense obstacles, and the target is a high-speed moving, highly flexible, strongly maneuverable and non-cooperative object, the limitations of traditional models will be greatly magnified, and even cause systematic failure.

[0003] Specifically, the complex operating environment and highly dynamic target pose serious challenges to the cooperative capture of UAV swarms. In an environment without a global positioning system, the UAV swarm cannot obtain the precise global coordinates of itself or its companions, and cannot accurately anchor the absolute positions of the target and obstacles. Therefore, the cooperative control strategy loses its spatial reference, which seriously restricts the situational awareness sharing and decision consistency of large-scale UAV swarms. In addition, the observation of a single UAV based on local sensors is easily affected by occlusion and noise interference, resulting in large deviations in estimation results and poor robustness. Traditional algorithms also significantly reduce accuracy when dealing with highly nonlinear motion, leading to misjudgments of the target's position, velocity and even intent. In order to establish global consistency between the UAV swarm and the target, simply fusing the observation data of multiple UAVs will face problems such as spatio-temporal misalignment and observation heterogeneity. Finally, in order to quickly enclose and capture the target, the UAV swarm must also meet the following conditions: quickly form an effective enclosure, strictly avoid obstacles, prevent intra-group collisions and maintain visibility to the target, etc., which are often conflicting constraints.

[0004] Traditional methods usually employ step-wise optimization or rule-based libraries, which makes it difficult to balance real-time performance and safety. This easily leads to a decision stalemate state of "oscillating wandering" or "path blocking" for the UAV. SUMMARY

[0005] In view of the above problems in the prior art, the unmanned aerial vehicle cluster cooperative target encirclement method based on brain-like computing provided by the present application solves the problem that the existing related methods are difficult to balance real-time performance and safety due to complex operation environment and highly dynamic target cooperative capture.

[0006] In order to achieve the above-mentioned purpose of the application, the technical scheme adopted by the present application is as follows: an unmanned aerial vehicle cluster cooperative target encirclement method based on brain-like computing, comprising the following steps:

[0007] S1, obtaining multi-view images under the camera view angle of the unmanned aerial vehicle, and obtaining the local pose and covariance of each unmanned aerial vehicle by using panoramic images, visual motion and grid-landmark cell model;

[0008] S2, obtaining the global consistent relative pose of the unmanned aerial vehicle cluster by using compression broadcast, co-view matching and graph optimization algorithm according to the local pose and covariance of each unmanned aerial vehicle;

[0009] S3, based on the global consistent relative pose of the unmanned aerial vehicle cluster, each unmanned aerial vehicle identifies the target, uses the unscented Kalman filter algorithm to estimate the state of the target, and fuses all target state estimates to obtain target state information in a unified coordinate system;

[0010] S4, based on the historical target state information, using the Bezier curve generation algorithm to obtain the target motion trajectory in a period of historical state, extending the Bezier curve trajectory to predict the motion of the target, and obtaining the predicted state of the target in the unified coordinate system;

[0011] S5, extending the unmanned aerial vehicle cluster by the method of centroid extension to obtain the encirclement queue for target encirclement, using the Gauss-Newton method to solve the minimum encirclement cost, and then obtaining the optimal encirclement queue to generate the encirclement position of each unmanned aerial vehicle;

[0012] S6, based on the encirclement position of each unmanned aerial vehicle, using the hybrid A* search algorithm to generate the motion trajectory of the unmanned aerial vehicle, performing anti-collision detection according to the motion trajectory, generating the final encirclement trajectory sequence, and completing the encirclement of the agile target.

[0013] Further, S1 comprises the following sub-steps:

[0014] S11, generating a continuous panoramic image by using the back projection and spherical mapping algorithm according to the multi-view images obtained by the monocular camera of the unmanned aerial vehicle cluster;

[0015] S12, based on the panoramic map, using ORB operator and RANSAC five-point method, the rotation matrix and translation vector between adjacent frames are obtained as the measured visual motion;

[0016] S13, based on the feature description amplitude, using the first pulse coding and grid-landmark cell model, position coding and memory are carried out;

[0017] S14, based on visual motion, IMU data and position coding, UKF is used to predict and correct the state, and the local pose and covariance of each UAV are output.

[0018] Further, in S11, the method of back projection and spherical mapping algorithm is: the multi-view image is back projected to the body coordinate system, wherein, for the first i UAV camera pixel The expression for back projection is:

[0019]

[0020]

[0021] In the formula, is the pixel coordinate back projected to the coordinate of the UAV body coordinate system, is the rotation matrix of the UAV in the world coordinate system, is the translation matrix of the starting point of the UAV in the world coordinate system, is the camera intrinsic matrix of the first i UAV, Y i is the coordinate of the body coordinate system, u is the horizontal coordinate of the image, v is the vertical coordinate of the image;

[0022] Map all Y i On the sphere with a radius to the angle space , and synthesize the panoramic map through bilinear interpolation ;

[0023] S12 is: according to the panoramic map, the key points and their descriptors are extracted using ORB operator, and are matched through FLANN, and the outliers are removed through RANSAC, and the relative motion is calculated under the known intrinsic parameters using five-point method, to obtain the measured visual motion ;

[0024]

[0025] In the formula, is the rotation matrix, It is a translation vector;

[0026] In S13, the grid-landmark cell model includes grid cells and landmark cell networks. Specifically, S13 maps feature description amplitudes to pulse times. , generating sparse pulses;

[0027]

[0028] In the formula, To prevent small amounts from being divided by zero, This is the upper limit of the encoding time window. For descriptors, The descriptor for the maximum amplitude at the current moment;

[0029] When sparse impulses are input into the grid cell and landmark cell network, the grid cell state vector is updated using the following formula;

[0030]

[0031] In the formula, It is a non-linear activation. and For learnable parameters, For time t The grid cell state vector, For time t +1 grid cell state vector;

[0032] Landmark cells activate associated memories by activating sparse pulse and grid cell state vectors according to the Hebbian rule;

[0033]

[0034] In the formula, To activate associated memories, To activate memories prior to association, For time t sparse pulses, For learning rate, It is the transpose symbol;

[0035] In subsequent frames, correction parameters are established based on the memory associated with activation. Correct the grid code, and use the corrected grid code as the position code;

[0036]

[0037] In the formula, For learning rate, For time t +1 sparse pulse;

[0038] In S14, the state vector of the UKF and observation are:

[0039]

[0040]

[0041] wherein, is the local pose of the time t , and is the velocity of the time t ;

[0042] The first estimate and the second estimate are obtained by the UKF nonlinear prediction, and the Kalman gain is calculated according to the second estimate;

[0043]

[0044] wherein, H is the state transition equation, is the noise matrix;

[0045] The state vector of the UKF is updated according to the Kalman gain and the first estimate, and the local pose and the covariance of each UAV are output by the UKF;

[0046]

[0047] wherein, is the updated state vector of the UKF, is the state transition equation transformation.

[0048] Further, S2 comprises the following sub-steps:

[0049] S21, based on the covariance and the sparse pulse, a broadcast data packet is generated for sharing the neighbor observation by using a differential quantization or an autoencoder algorithm;

[0050] S22, the field of view overlap area is estimated according to the broadcast data packet, and the relative observation between the UAVs is calculated by using an ORB+RANSAC algorithm;

[0051] S23, based on the relative observation between the UAVs, the global pose and the covariance of the UAVs are solved by using an incremental Gauss-Newton algorithm;

[0052] S24, the global pose and the covariance of the UAVs are sent to each UAV, the global pose is corrected in a closed loop, and the global consistent relative pose of the UAV cluster is obtained.

[0053] Further, S21 is specifically: the local pose increment and covariance output by each UAV are differentially quantized to obtain differentially quantized parameters and :

[0054]

[0055]

[0056] wherein, is a unified quantization function of vectors and matrices, is the covariance of time t -1, is the covariance of time t +1, is the local pose increment;

[0057]

[0058] wherein, is the local pose of time t -1,

[0059] The grid cell state vector and sparse pulses are dimensionally reduced by using a sparse autoencoder to obtain a dimensionally reduced result h ;

[0060]

[0061]

[0062] wherein, is a spliced vector, is an activation function, is a sparse autoencoder parameter, is a state vector corresponding to an estimated value;

[0063] A broadcast data packet is established , and the broadcast data packet reaches the adjacent UAV within a maximum delay allowed by the link quality;

[0064] S22 is specifically: the adjacent UAV grid activation is obtained by the UAV receiving the broadcast data packet of the adjacent UAV, the visual field center deviation is calculated according to the UAV self grid activation and the adjacent UAV grid activation, the visual field overlapping center and its radius are estimated, and the visual field overlapping area is obtained;

[0065] The visual field center deviation includes a first spherical coordinate deviation and a second spherical coordinate deviation ;

[0066]

[0067]

[0068] In the formula, Activate the drone's own mesh. Activate the neighboring machine's mesh. and This is a mapping function from grid encoding to spherical coordinate differences;

[0069] ORB operators are used to extract segments of the decoded panoramic image in the overlapping regions of the field of view, and joint descriptors are constructed by combining the temporal information of the first pulse encoding. Matching pairs are obtained after RANSAC filtering. In the known camera intrinsic parameters K Solving the relative observations between UAVs under certain conditions This includes the rotation matrix of the drone itself relative to its neighboring drones. Translation vector ;

[0070]

[0071] In the formula, R For the first i The drone and the first j The difference in rotational variation between the drones For the first i A drone in the camera's internal reference K In this case, the difference in translation from the initial position, For the first j A drone in the camera's internal reference K In this case, the difference in translation from the initial position;

[0072] Observation covariance of relative observations between UAVs The specific expression is:

[0073]

[0074] In the formula, The number of all drones to establish contact;

[0075] S23 specifically involves: establishing a graph optimization model for a drone swarm, where the node set represents the local pose of each drone, the edge set consists of the relative observations and observation covariance among the drones, and the objective function of the graph optimization model is... The specific expression is:

[0076]

[0077] In the formula, For the first in three-dimensional space jThe difference in pose between the drone and its initial position. For the first in three-dimensional space i The difference in pose between the drone and its initial position. for and The pose difference between them For SE(3) complex;

[0078] The incremental Gauss-Newton algorithm works as follows: upon receiving new relative observations, it linearizes the residuals. The sparse information matrix is ​​incrementally updated, and the local pose of the UAV is updated after the increment is solved, which serves as the global pose of the UAV.

[0079]

[0080] In the formula, For the objective function The first partial derivative, For the objective function The first differential, For the objective function The second partial derivative, For the objective function The second differential, This is the difference between the sum of the partial derivative and the differential product;

[0081]

[0082] In the formula, The sparse information matrix before the update. This is the updated sparse information matrix;

[0083]

[0084] In the formula, p i This is the local pose of the drone before the update. This represents the updated local pose of the drone.

[0085] S24 specifically involves broadcasting the global pose and covariance of the UAV to each UAV, and updating the UKF state of the local UAV as follows:

[0086]

[0087]

[0088] In the formula, The updated state vector, This represents the global pose of the drone. This is the updated grid cell state vector.P i For covariance, and For the prior noise of velocity and mesh activation, To update the pose of the drone in three-dimensional space. For the updated speed;

[0089] Various drones utilize and Restart the filtering loop and input the updated grid cell state vector and sparse impulses into the landmark cells for long-term memory updates to obtain the globally consistent relative pose of the UAV cluster.

[0090] Furthermore, S3 includes the following sub-steps:

[0091] S31. Using the YOLO target recognition method, identify the target in the image information and obtain the target position in its coordinate system based on the image information;

[0092] S32. Transform the target position using a rotation matrix to obtain the target position in a unified coordinate system;

[0093] S33. Based on the target position in a unified coordinate system, UKF is used to estimate the target state, and the target state estimate and covariance after filtering out sensor errors are obtained.

[0094] S34. Using a loosely coupled approach, the noise covariance of each UAV's observation of the target is combined to fuse all target state estimates and obtain target state information.

[0095] S33 specifically refers to:

[0096] Based on the state vector using UKF Covariance Matrix Generate 2 n The formula for generating +1 sigma point is as follows:

[0097]

[0098] In the formula, For sigma points, n Let be the dimension of the state vector. for K The estimated value corresponding to the position at a given time. It is a regulation function;

[0099]

[0100] In the formula, and To adjust the parameters;

[0101] Propagate sigma points through the nonlinear state transition function to compute the predicted state mean and covariance ;

[0102]

[0103]

[0104] where is the predicted value of the state at time K +1 from the state value at time K +1, Q is the covariance matrix of the measurement noise, , is the state part, is the noise part, and are weight coefficients;

[0105]

[0106]

[0107] where is the small error;

[0108] Map the predicted Sigma points to the measurement space to compute the mean , covariance matrix S and cross-covariance matrix T ;

[0109]

[0110]

[0111]

[0112] where is the predicted value of the observation at time K +1 from the observation value at time K +1, R i is the assumed predetermined noise covariance matrix;

[0113]

[0114] Compute the Kalman gain from the covariance matrix ;

[0115]

[0116] Compute the innovation from the observation and the measured mean value updating the state estimation and covariance of the target, the updated target state estimation and covariance The expression is:

[0117]

[0118]

[0119] In the formula, is K observation value at the moment;

[0120] In S34, the target state information is obtained The expression is specifically:

[0121]

[0122] In the formula, is the target state estimation of the i unmanned aerial vehicle, is the weight of the unmanned aerial vehicle;

[0123]

[0124] In the formula, is the observation noise covariance of the unmanned aerial vehicle.

[0125] Further, S4 includes the following sub-steps:

[0126] S41, using a sliding window method to record historical target state information of a set length;

[0127] S42, based on the historical target state information, generating a target motion trajectory in the historical state according to the Bezier curve generation algorithm;

[0128] S43, calculating the predicted step length by the hyperbolic tangent function according to the distance between the center of mass of the unmanned aerial vehicle cluster and the current target position;

[0129] S44, extending the predicted step length of the Bezier curve to obtain the predicted state of the target in the unified coordinates;

[0130] S42 specifically: using a Bezier curve to describe the target motion trajectory, and the target motion trajectory The expression is specifically:

[0131]

[0132] In the formula, is an n-th Bernstein polynomial base, is a control point of the Bezier curve;

[0133] In S43, the predicted step length is calculated The expression is specifically:

[0134]

[0135] In the formula, is a preset predicted time step, d is a target current position distance, is a current time.

[0136] Further, S5 includes the following sub-steps:

[0137] S51, on the basis of the global unified coordinates, the centroid of the UAV group is calculated;

[0138] S52, the vector of each UAV from the centroid to the UAV is calculated, all UAVs are expanded to form a surrounding capture queue, the UAVs in the centroid area and the small sector area are sorted in a permutation and combination manner to form a plurality of surrounding capture queues;

[0139] S53, the Gauss-Newton method is used to solve the surrounding capture function to obtain the minimum surrounding capture cost corresponding to the surrounding capture queue;

[0140] S54, the minimum surrounding capture costs corresponding to all candidate surrounding capture queues are compared, the surrounding capture queue with the smallest minimum surrounding capture cost is taken as the optimal surrounding capture queue, the optimal surrounding capture point is calculated to obtain the surrounding capture position of each UAV;

[0141] S52 is specifically: the vector of each UAV relative to the centroid is obtained, and the UAV included angle of the vector and the surrounding capture direction vector is calculated, wherein the expression for calculating the included angle i of the first UAV is specifically:

[0142]

[0143] In the formula, is a modulo function, is the included angle of the vector and the axis, and the value range is [0, 2π), x is the vector of the first UAV relative to the centroid, i , is the centroid, C , is the number of UAVs, N is the surrounding capture direction vector, , is the target surrounding capture point; T ​

[0144] Sort all the angles of the UAVs, and assign a serial number to each UAV according to the sorting result to generate a preliminary encirclement sequence;

[0145] If there is a UAV at the centroid position, sequentially insert the UAV between the two adjacent UAVs in the circular formation, and generate an encirclement queue according to all possible insertion sequences;

[0146] If there are more than a set number of UAVs in the same direction, and the change of the angles of the UAVs in the same sector is less than a preset threshold, then arrange and combine the set of UAVs to obtain all possible insertion sequences, and insert them into the interval positions of the preliminary encirclement sequence to generate an encirclement queue according to all possible insertion sequences;

[0147] In S53, the minimum encirclement cost corresponding to the encirclement queue is calculated The expression is specifically:

[0148]

[0149] In the formula, S is the encirclement queue, is the initial angular offset, R is the encirclement radius, is the cost function;

[0150]

[0151] In the formula, is the actual position of the UAV in the encirclement queue S . j

[0152] Further, S6 includes the following sub-steps:

[0153] S61, obtain surrounding obstacle information from the binocular camera, construct a three-dimensional grid or octree environment model, and map the continuous state of the UAV to a discrete grid and a heading discrete angle, which is used for searching on the grid map. The continuous state includes position, orientation, and velocity;

[0154] S62, take the current position, orientation, and velocity of the UAV as the search starting point, and determine the safe encirclement point and posture of the target region or target swarm according to the predetermined "capture encirclement" strategy;

[0155] S63, in the discretized state space, take the kinematic constraints of the UAV as the expansion rule, use a heuristic function to evaluate the cost of the current grid point to the target, and promote the search to develop towards the target direction, which is used to generate a trajectory that satisfies the dynamics constraints and is connected;

[0156] ​S64, spline fitting or smoothing processing of the discrete trajectory generated based on the hybrid state A algorithm to eliminate sharp turns of the polyline segment, assigning speed and time stamp to the smooth trajectory according to the speed and acceleration limit of the unmanned aerial vehicle, to obtain a continuous executable time parameterized trajectory;

[0157] S65, discretely sampling the smooth trajectory, using the nearest obstacle distance detection algorithm to detect obstacle collision at each sampling point, and in response to detecting collision or insufficient safety distance, triggering local re-planning to perform high-low flight near the segment to correct the path in the smooth trajectory;

[0158] S66, integrating and splicing the smooth trajectory, the time parameterized trajectory and the trajectory segment that has passed the collision detection, to generate a final pursuit trajectory sequence;

[0159] S64 specifically is:

[0160] Target-oriented kinematic search front end: constructing a search tree based on the hybrid state A algorithm, generating motion primitives by discretizing control input, and introducing a cost function in the search process and the target motion trajectory B(t) as heuristic information, designing a double heuristic function:

[0161]

[0162] In the formula, is a dynamic distance estimation based on the optimal boundary value problem, is the target state of the unmanned aerial vehicle, is the current state of the unmanned aerial vehicle, is the weight, S t is the sum of the expected expansion time, is the path expansion time, obtained by extrapolating the predicted trajectory:

[0163]

[0164] In the formula, is the Bezier curve trajectory of the original historical trajectory, is the trajectory obtained by time extrapolation;

[0165] Space-time trajectory optimization back end: based on the initial path generated by the target-oriented kinematic search front end, constructing a flight corridor composed of connected free space cubes ;

[0166] In the formula, is the Bezier curve control point, M is the number of trajectory segments of the flight corridor;

[0167] Establish an optimization problem with dynamic constraints under corridor constraints:

[0168]

[0169] In the formula, For smoothing terms, For the field of power, For dynamic constraint terms, q w intermediate waypoints T For the time of segmented trajectory, for x , y and z The combination x For segmented waypoints x Axis coordinates y For segmented waypoints y Axis coordinates z For segmented waypoints z Axis coordinates; This optimization problem is solved efficiently using the split convex optimization method, generating a time-parameterized trajectory;

[0170] In S66, the collision detection method is as follows:

[0171] For each UAV's time-parameterized trajectory, let the UAV... The trajectory is m Composition of segmental polynomials, defining the parameter set ;

[0172]

[0173] In the formula, For the first k Duration of segment , For the real number field, For the first k Segment polynomial order , For the matrix field, The coefficient matrix, , For the first k segment polynomial x Parameters in the axial direction, For the first k segment polynomial y Parameters in the axial direction, For the first k segment polynomial z Parameters in the axial direction;

[0174] Calculate the relative time of the current moment within each segment using a polynomial. ;

[0175]

[0176] wherein, is the difference between the current time and the start time, is the absolute time, is the track start time stamp, is the maximum function, s is the number of track segments, is the minimum index function, is the time of each track segment;

[0177] calculating the three-dimensional space position of the predicted track ;

[0178]

[0179] wherein, is the current position of the UAV, is the coefficient of the s th segment, j th item, , is the coefficient in the x axis direction, is the coefficient in the y axis direction, is the coefficient in the z axis direction, is the relative time within the s th segment;

[0180] comparing the sampling point of each time in the predicted track of the UAV with the position of the predicted track of other UAVs at the same time, if the distance is lower than the preset safety distance, it is determined that there is a potential collision risk;

[0181] when it is determined that there is a potential collision risk, each UAV with a potential collision risk receives the state information of the target in real time, and obtains the position information of other UAVs, each UAV calculates the Euclidean distance between itself and the target, and compares the Euclidean distance with the Euclidean distance between the target and other UAVs, sorts the UAVs in ascending order according to the Euclidean distance to determine the motion priority, the UAV with the highest priority keeps the original planned track, and the other UAVs avoid obstacles at different heights according to the priority, and the expression of the sequence number after the priority sorting is:

[0182]

[0183] wherein, is the position coordinate of the j th UAV, Target position coordinates, rank(•) is a ranking function.

[0184] The beneficial effects of the present application are:

[0185] (1) The method solves the problem of relative positioning of the UAV group in the rejection environment, and the global positioning is constructed through the cluster cooperative fusion and global optimization, thereby laying a foundation for cooperative operation, and the UAV group can effectively cope with environmental uncertainty and its own maneuverability.

[0186] (2) The method makes changes in the cooperative target state estimation and prediction algorithm, and fuses the unscented Kalman filter algorithm and the Bezier curve generation algorithm in the multi-source observation trajectory prediction method, which significantly improves the accuracy of motion state estimation.

[0187] (3) The method makes changes in the hunting strategy algorithm, and innovatively proposes a centroid extension optimization strategy generation algorithm based on virtual centroid, as well as a space-time trajectory optimization method, to realize fast target surrounding under the premise of ensuring safety and generate a satisfactory dynamic feasible trajectory. BRIEF DESCRIPTION OF DRAWINGS

[0188] Figure 1 A flow chart of a UAV cluster cooperative target hunting method based on brain-like computing. DETAILED DESCRIPTION

[0189] The specific embodiments of the present application are described below to facilitate understanding of the present application by those skilled in the art, but it should be clear that the present application is not limited to the scope of the specific embodiments. For those skilled in the art, it is obvious that various changes are within the spirit and scope of the present application defined and determined by the appended claims, and all inventions utilizing the concept of the present application are within the scope of protection.

[0190] As shown in Figure 1 In one embodiment of the present application, a UAV cluster cooperative target hunting method based on brain-like computing includes the following steps:

[0191] S1, obtaining multi-view images under camera view angles according to UAVs, using panoramic images, visual motion and grid-landmark cell model to obtain local pose and covariance of each UAV;

[0192] S2, obtaining global consistent relative pose of the UAV cluster through compression broadcast, co-view matching and graph optimization algorithm according to the local pose and covariance of each UAV;

[0193] S3. Based on the globally consistent relative pose of the UAV swarm, each UAV identifies the target, uses the unscented Kalman filter algorithm to estimate the target state, and fuses all the target state estimates to obtain the target state information in a unified coordinate system.

[0194] S4. Based on historical target state information, use the Bézier curve generation algorithm to obtain a target motion trajectory under a historical state, expand the Bézier curve trajectory to predict the target motion, and obtain the target prediction state under a unified coordinate system.

[0195] S5. Expand the drone swarm by extending the centroid to obtain the target capture queue. Use the Gauss-Newton method to solve for the minimum capture cost, and then obtain the optimal capture queue to generate the capture position of each drone.

[0196] S6. Based on the capture position of each drone, a hybrid A* search algorithm is used to generate the drone's motion trajectory. Collision avoidance detection is performed based on the motion trajectory to generate the final capture trajectory sequence, thus completing the capture of the agile target.

[0197] S1 includes the following steps:

[0198] S11. Based on the multi-view images obtained by the UAV cluster through the monocular camera, a continuous panoramic image is generated using the back projection and spherical mapping algorithm.

[0199] S12. Based on the panoramic image, the ORB operator and RANSAC five-point method are used to obtain the rotation matrix and translation vector between adjacent frames, which are used as the measurement of visual motion.

[0200] S13. Based on the feature description amplitude, the first pulse coding and grid-landmark cell model are used to perform position coding and memory;

[0201] S14. Based on visual motion, IMU data and position encoding, the UKF (Unscented Kalman Filter) algorithm is used to predict and correct the state, and output the local pose and covariance of each UAV.

[0202] In S11, the method using back projection and spherical mapping algorithm is as follows: the multi-view image is back-projected onto the body coordinate system, wherein, for the ... i drone camera pixels The specific expression for back projection is:

[0203]

[0204]

[0205] In the formula, These are the pixel coordinates projected onto the UAV's body coordinate system. is a rotation matrix of the UAV in the world coordinate system, is a translation matrix of the UAV in the world coordinate system at the starting point (origin), is the first i is the camera intrinsic matrix of the UAV, Y i is the coordinate of the body coordinate system, u is the coordinate of the image in the horizontal direction, v is the coordinate of the image in the vertical direction;

[0206] All Y i is mapped to the angle space on a spherical surface with a radius , and a panoramic image is synthesized by bilinear interpolation ; in this embodiment, the process of synthesizing the panoramic image maintains the geometric consistency of each camera view, and provides all-around environmental observation for subsequent visual processing.

[0207] S12 is specifically: extracting key points and their descriptors using the ORB operator according to the panoramic image, and pairing through FLANN, removing outliers by RANSAC, and calculating the relative motion under the known intrinsic parameter by the five-point method to obtain the measured visual motion ;

[0208]

[0209] In the formula, is a rotation matrix, is a translation vector;

[0210] In S13, the grid-landmark cell model includes a grid cell and a landmark cell network, and S13 is specifically: mapping the feature description amplitude to the pulse time to generate sparse pulses;

[0211]

[0212] In the formula, is a small amount to prevent division by zero, is the upper limit of the encoding time window, is the descriptor, is the descriptor of the maximum amplitude at the current time;

[0213] The sparse pulses are input into the grid cell and the landmark cell network, and the grid cell state vector is updated by the following formula;

[0214]

[0215] In the formula, is a nonlinear activation, and is a learnable parameter, is a time t grid cell state vector, is a time t +1 grid cell state vector;

[0216] The landmark cell activates the sparse pulse and the grid cell state vector to associate the memory according to the Hebbian rule;

[0217]

[0218] wherein, is the activated associated memory, is the memory before the activation association, is a time t sparse pulse, is a learning rate;

[0219] The correction parameter is established according to the activated associated memory in the subsequent frame corrects the grid code, and the corrected grid code is used as the position code;

[0220]

[0221] wherein, is a transpose symbol, is a learning rate, is a time t +1 sparse pulse;

[0222] In S14, the state vector and the observation of the UKF are:

[0223]

[0224]

[0225] wherein, is a time t local pose, is a time t velocity;

[0226] The first estimation and the second estimation are obtained through the UKF nonlinear prediction, and the Kalman gain is calculated according to the second estimation;

[0227]

[0228] wherein, H is a state transition equation, This is the noise matrix;

[0229] The UKF state vector is updated based on the Kalman gain and the first prediction, and the local pose and covariance of each UAV are output through the UKF.

[0230]

[0231] In the formula, This is the UKF-updated state vector. This is a transformation of the state transition equation.

[0232] In this embodiment, the final output is the local pose. With covariance This provides a reliable input for cluster integration.

[0233] S2 includes the following steps:

[0234] S21. Based on covariance and sparse impulses, differential quantization or autoencoder algorithms are used to generate broadcast data packets for sharing neighboring machine observations.

[0235] S22. Estimate the overlapping area of ​​the field of view based on the broadcast data packets, and use the ORB+RANSAC algorithm to calculate the relative observations between UAVs;

[0236] S23. Based on the relative observation between UAVs, the incremental Gauss-Newton algorithm is used to solve the global pose and covariance of the UAVs, and the global pose is solved in real time.

[0237] S24. The global pose and covariance of the UAV are sent to each UAV, and the global pose is corrected in a closed loop to obtain the globally consistent relative pose of the UAV cluster, thereby improving the closed-loop accuracy.

[0238] S21 specifically involves: performing differential quantization on the output local pose increment and covariance of each UAV to obtain the differentially quantized parameters. and :

[0239]

[0240]

[0241] In the formula, A unified quantization function for vectors and matrices. For time t covariance, For time t +1 covariance, For local pose increments;

[0242]

[0243] wherein, is the time t at which the local pose is

[0244] In this embodiment, to reduce bandwidth overhead, each UAV differentially quantizes the local pose increment and the covariance of its output.

[0245] The grid cell state vector and sparse pulses are dimensionally reduced using a sparse autoencoder to obtain the reduced dimension result h .

[0246]

[0247]

[0248] wherein, is the stitching vector, is an activation function, is a sparse autoencoder parameter, is a state vector corresponding to the estimated value;

[0249] A broadcast data packet is established , and the broadcast data packet reaches the adjacent UAV within a maximum delay allowed by the link quality;

[0250] S22 is specifically: the UAV receives the broadcast data packet of the adjacent UAV to obtain the grid activation of the adjacent UAV, calculates the visual center deviation according to the grid activation of the UAV itself and the grid activation of the adjacent UAV, estimates the visual overlapping circle center and its radius, and obtains the visual overlapping area;

[0251] The visual center deviation includes a first spherical coordinate deviation and a second spherical coordinate deviation .

[0252]

[0253]

[0254] wherein, is the grid activation of the UAV itself, is the grid activation of the adjacent UAV, and is a mapping function from the grid encoding to the spherical coordinate difference value, which is used to estimate the visual overlapping circle center and its radius, so as to obtain the visual overlapping area;

[0255] ORB operators are used to extract segments of the decoded panoramic image in the overlapping regions of the field of view, and joint descriptors are constructed by combining the temporal information of the first pulse encoding. Matching pairs are obtained after RANSAC filtering. In the known camera intrinsic parameters K Solving the relative observations between UAVs under certain conditions This includes the rotation matrix of the drone itself relative to its neighboring drones. Translation vector ;

[0256]

[0257] In the formula, R For the first i The drone and the first j The difference in rotational variation between the drones For the first i A drone in the camera's internal reference K In this case, the difference in translation from the initial position, For the first j A drone in the camera's internal reference K In this case, the difference in translation from the initial position;

[0258] Observation covariance of relative observations between UAVs The specific expression is:

[0259]

[0260] In the formula, The number of all drones to establish contact;

[0261] S23 specifically involves: establishing a graph optimization model for a drone swarm, where the node set represents the local pose of each drone, the edge set consists of the relative observations and observation covariances among the drones, and the objective function of the graph optimization model is... The specific expression is:

[0262]

[0263] In the formula, For the first in three-dimensional space j The difference in pose between the drone and its initial position. For the first in three-dimensional space i The difference in pose between the drone and its initial position. for and The pose difference between them For SE(3) complex;

[0264] The working process of the incremental Gauss-Newton algorithm is as follows: when receiving new relative observations, linearize the residual , and incrementally update the sparse information matrix to solve the increment , and update the local pose of the UAV after the update as the global pose of the UAV;

[0265]

[0266] In the formula, is the first partial derivative of the target function , is the first differential of the target function , is the second partial derivative of the target function , is the second differential of the target function , is the difference between the sum of the product of the partial derivative and the differential;

[0267]

[0268] In the formula, is the sparse information matrix before the update, is the sparse information matrix after the update;

[0269]

[0270] In the formula, p i is the local pose of the UAV before the update, is the local pose of the UAV after the update;

[0271] S24 is specifically: broadcasting the global pose and covariance of the UAV to each UAV, and the UKF state update of the local UAV is specifically:

[0272]

[0273]

[0274] In the formula, is the updated state vector, is the global pose of the UAV, is the updated grid cell state vector, P i is the covariance, and are the prior noise of the velocity and the grid activation, is the updated pose in the three-dimensional space of the UAV, is the updated velocity;

[0275] ​​​​Each UAV utilizes and restarts the filtering cycle, and updates the long-term memory of the updated grid cell state vector and sparse pulse input landmark cell, to obtain the global consistent relative pose of the UAV cluster.

[0276] Through the above closed-loop correction, the system tightly couples the local and global solutions in large-scale formation, and realizes high-precision and low-drift end-to-end relative positioning.

[0277] S3 includes the following sub-steps:

[0278] S31, a YOLO target recognition method is used to recognize the target in the image information, and the target position in the coordinate system thereof is obtained according to the image information;

[0279] S32, the target position is transformed through a rotation matrix to obtain the target position in the unified coordinate system;

[0280] S33, based on the target position in the unified coordinate system, UKF is used to estimate the state of the target to obtain the target state estimation and covariance after filtering the sensor error;

[0281] S34, in a loosely coupled manner, the observation noise covariance of each UAV on the target is combined to fuse all target state estimations to obtain the target state information;

[0282] S33 specifically includes:

[0283] 2 +1 sigma points are generated by UKF according to the state vector and the covariance matrix n , and the generation formula is specifically as follows:

[0284]

[0285] In the formula, is a sigma point, n is the dimension of the state vector, is the estimated value corresponding to the position at the moment, K is an adjustment function;

[0286]

[0287] In the formula, and are adjustment parameters;

[0288] The sigma points are propagated through a nonlinear state transition function to calculate the predicted state mean and covariance ;​

[0289]

[0290]

[0291] wherein, is predicted by K the state value at time K +1, Q is a covariance matrix of measurement noise, , is a state part, is a noise part, and is a weight coefficient;

[0292]

[0293]

[0294] wherein, is a small error;

[0295] mapping the predicted Sigma points to the measurement space, calculating the mean value of the measurement , the covariance matrix S and the cross-covariance matrix T , calculating the cross-covariance matrix for determining the accuracy of the calculated measurement value;

[0296]

[0297]

[0298]

[0299] wherein, is predicted by K the observation value at time K +1, R i is a predetermined noise covariance matrix assumed;

[0300]

[0301] calculating the Kalman gain from the covariance matrix;

[0302]

[0303] updating the state estimation and the covariance of the target from the observation value and the mean value of the measurement the updated target state estimation and covariance The expression is:

[0304]

[0305]

[0306] In the formula, is the observation value at the moment; K

[0307] In S34, the expression of the target state information is specifically:

[0308]

[0309] In the formula, is the target state estimation of the jth unmanned aerial vehicle, i is the weight of the jth unmanned aerial vehicle;

[0310]

[0311] In the formula, is the observation noise covariance of the jth unmanned aerial vehicle. In the embodiment, multiple unmanned aerial vehicles observe the position of the target, and a weighted average method is used to fuse the observation results of each unmanned aerial vehicle to obtain the target state information. Through the method, the fast sharing of the accurate target state information of the cluster can be realized.

[0312] S4 includes the following sub-steps:

[0313] S41, a sliding window method is used to record historical target state information of a set length;

[0314] S42, based on the historical target state information, a target motion trajectory in the historical state is generated according to the Bezier curve generation algorithm;

[0315] S43, according to the distance between the center of mass of the unmanned aerial vehicle cluster and the current target position, the predicted step is calculated through the hyperbolic tangent function;

[0316] S44, the Bezier curve is extended by the predicted step to obtain the target predicted state in the unified coordinate;

[0317] S42 is specifically: the Bezier curve is used to describe the target motion trajectory, and the expression of the target motion trajectory

[0318]

[0319] ​​​​​

[0320] wherein, is an n-th Bernstein polynomial basis, is a Bezier curve control point;

[0321] In the embodiment, the target motion trajectory is described by Bernstein basis polynomials, referred to as B´ezier curves, to describe the target prediction trajectory. The present application takes t∈ The 3D position of the target observed in the global frame at time t is denoted as Then a FIFO queue with length L is maintained to store past observations and corresponding timestamps. The queue is denoted as [ , ,…, ] where { , }. The time range contained is [ , ] where is equal to the current time. When a new target observation is obtained, a new target motion trajectory is generated by fitting past observations.

[0322] In S43, the expression of the predicted step is calculated as follows:

[0323]

[0324] wherein, is a preset prediction time step, d is the distance between the target current position and the centroid, is the current time.

[0325] In the embodiment, the present application uses the hyperbolic tangent function tanh(x) to determine the prediction target position of the pursuit when the distance d between the centroid and the target current position is different, in order to quickly pursue the unintended target during the pursuit.

[0326] S5 includes the following steps:

[0327] S51, calculate the centroid of the UAV swarm on the basis of the global unified coordinates;

[0328] S52, calculate the vector from the centroid to each UAV, and all UAVs expand outward to form a pursuit queue. The UAVs in the centroid area and the small sector area are sorted in a permutation and combination manner to form several pursuit queues.

[0329] S53, solve the trapping function by using Gauss-Newton method to obtain the minimum trapping cost corresponding to the trapping queue;

[0330] S54, compare the minimum trapping costs corresponding to all candidate trapping queues, take the trapping queue with the minimum minimum trapping cost as the optimal trapping queue, calculate the optimal trapping point, and thus obtain the trapping position of each UAV;

[0331] In S51, at time t, the position of the UAV group is known, the position of each UAV is , , N , T , T , ,

[0332] S52 is specifically: obtain the vector of each UAV relative to the centroid, calculate the UAV included angle between the vector and the trapping direction vector, wherein the expression for calculating the included angle of the first i UAV is ,

[0333]

[0334] In the formula, , is the included angle between the vector and the axis, and the value range is [0, 2π), x , is the vector of the first i UAV relative to the centroid, , C is the centroid;

[0335] Sort all UAV included angles, assign serial numbers to UAVs according to the sorting result, and generate a preliminary trapping order;

[0336] In this embodiment, ,wherein θ( ) represents the included angle between the vector and the axis, and the value range is [0, 2π); here, the sorting is in the clockwise direction. Sort all x , and reassign serial numbers to the sorted UAVs, satisfying .

[0337] If there is a UAV in the centroid position, the UAV is sequentially inserted between two adjacent UAVs in the circular formation, and a surrounding formation queue is generated according to all possible insertion arrangement orders. Based on this, all possible insertion arrangement orders Π(J) are obtained, and all candidate sequences thus obtained are written into a FIFO candidate queue

[0338] If there are more than a set number of UAVs in the same direction, and the change of the included angle of the dense UAVs in the same sector is less than a preset threshold, the UAV set is arranged and combined to obtain all possible insertion orders, which are inserted into the interval position of the preliminary surrounding order, and a surrounding formation queue is generated according to all possible insertion arrangement orders. Based on this, all possible insertion arrangement orders Π(J) are obtained, and a certain arrangement σ(J) is selected from them, which is inserted into the interval position of the original sequence to form a candidate adjacent sequence, and all candidate sequences thus obtained are written into a FIFO candidate queue

[0339] In S53, the expression of the minimum surrounding cost corresponding to the surrounding formation queue is calculated The expression is specifically as follows:

[0340]

[0341] In the formula, S is the surrounding formation queue, is the initial angular offset, R is the surrounding radius, is the cost function; and the minimum surrounding cost The derivation steps of the minimum surrounding cost are specifically as follows:

[0342] The expected position of each UAV on the ideal circular surrounding formation is set as:

[0343]

[0344] In the formula, R is the surrounding radius, is the initial angular offset, which is an optimization parameter.

[0345] The surrounding cost is calculated for each adjacent case (i.e., a sorting sequence) in the candidate queue The cost function is defined as:

[0346]

[0347] In the formula, is the surrounding formation queue S , the i-th j ​​actual position of the drone. Next, the parameters are solved for each candidate adjacency case by the Gauss-Newton method. By using the second-order expansion of Taylor series, the minimum value of the function is solved by iteration through the gradient and Hessian matrix. For the objective function f(x), the iteration formula of the Newton method is as follows:

[0348]

[0349] where, is the parameter estimation value at the k th iteration, is the Hessian matrix, and the cost function to be minimized is The cost function is taken with respect to and R, respectively, and the gradient is obtained as follows:

[0350]

[0351]

[0352] The elements of the Hessian matrix are the second-order partial derivatives of the cost function, which can be expressed as:

[0353]

[0354] The parameters are updated using the Gauss-Newton method iteration formula, and the specific expression is as follows:

[0355]

[0356] If the modulus of the gradient is less than a given threshold, it is considered to be converged, and the iteration is stopped. The corresponding minimum trapping cost is calculated as:

[0357]

[0358] S6 includes the following sub-steps:

[0359] S61, obtain the surrounding obstacle information from the binocular camera, construct a three-dimensional grid or octree environment model, map the continuous state of the drone to discrete grid points and heading discrete angles, for searching on the grid point graph, the continuous state includes position, orientation and velocity;

[0360] S62, take the current position, orientation and velocity of the drone as the search starting point, determine the safe trapping point and attitude of the target region or target swarm according to the predetermined "trapping" strategy;

[0361] ​S63, in the discretized state space, taking the kinematic constraints of the UAV as extension rules, using heuristic functions to evaluate the cost of the current grid point to the target, and pushing the search towards the target direction, for generating trajectories that satisfy the kinematic constraints and are connected;

[0362] S64, performing spline fitting or smoothing processing of the minimum bending energy on the discrete trajectory generated based on the hybrid state A algorithm, eliminating the sharp turns of the polyline segment, assigning velocity and time stamp to the smoothed trajectory according to the velocity and acceleration limits of the UAV, and obtaining a continuous executable time parameterized trajectory;

[0363] S65, discretely sampling the smoothed trajectory, using the nearest obstacle distance detection algorithm to detect obstacle collision at each sampling point, and in response to detecting collision or insufficient safety distance, triggering local re-planning to execute high-low air in the vicinity of the segment to correct the path in the smoothed trajectory;

[0364] S66, integrating and splicing the smoothed trajectory, the time parameterized trajectory and the trajectory segment that has passed the collision detection, to generate a final pursuit trajectory sequence;

[0365] S64 specifically includes:

[0366] Target-oriented kinematic search front end: constructing a search tree based on the hybrid state A algorithm, generating motion primitives by discretizing control inputs, and introducing a cost function in the search process and the target motion trajectory B(t) as heuristic information, designing a double heuristic function:

[0367]

[0368] wherein, is the dynamic distance estimation based on the optimal boundary value problem, is the target state of the UAV, is the current state of the UAV, is the weight, S t is the sum of the expected expansion time, is the path expansion time, The prediction trajectory is obtained by extrapolation:

[0369]

[0370] wherein, is the Bezier curve trajectory of the original historical trajectory, is the trajectory obtained by time extrapolation; this design makes the search process have foresight of the target motion, significantly improving the search efficiency in complex environments.

[0371] Space-time trajectory optimization backend: build a flight corridor consisting of connected free-space cubes based on the initial path generated by the goal-oriented kinematic search frontend ;

[0372] wherein, is a Bezier curve control point, M is the number of trajectory segments of the flight corridor;

[0373] An optimization problem with dynamics constraints is established under the corridor constraint:

[0374]

[0375] wherein, is a smoothing term, is a potential field term, is a dynamic constraint term, q w is an intermediate waypoint, T is the time of the segmented trajectory, is x , y and z a combination of, x is the x axis coordinate in the segmented waypoint, y is the y axis coordinate in the segmented waypoint, z is the z axis coordinate in the segmented waypoint; the optimization problem is efficiently solved by a split convex optimization method to generate a time-parameterized trajectory;

[0376] In this embodiment, the smoothing term adopts third derivative integral minimization, the potential field term ensures trajectory safety through a logarithmic obstacle function, and the dynamic constraint term adopts a segmented linear penalty function to process speed and acceleration constraints. The optimization problem is efficiently solved by a split convex optimization method to generate a dynamically feasible trajectory that satisfies.

[0377] In S66, under high-speed pursuit work, collision detection is performed on the next moment of the predicted path of the UAV to reduce the amount of calculation and improve real-time performance. The collision detection method is specifically:

[0378] For the time-parameterized trajectory of each UAV, the trajectory of the UAV is composed of m segmented polynomials, and the parameter set is defined.

[0379]

[0380] wherein, is the duration of the k segment, , For the real number field, For the first k Segment polynomial order , For the matrix field, The coefficient matrix, , For the first k segment polynomial x Parameters in the axial direction, For the first k segment polynomial y Parameters in the axial direction, For the first k segment polynomial z Parameters in the axial direction;

[0381] Calculate the relative time of the current moment within each segment using a polynomial. ;

[0382]

[0383] In the formula, It is the difference between the current time and the starting time. For absolute time, This is the starting timestamp of the trajectory. It is a function with maximum value. s The number of trajectory segments, For the minimum value index function, The time for each trajectory segment;

[0384] Calculate the three-dimensional spatial position of the predicted trajectory ;

[0385]

[0386] In the formula, This is the current location of the drone. For the first s Duan Di j Term coefficient, , for x Coefficient in the axial direction, for y Coefficient in the axial direction, for z Coefficient in the axial direction, In the first s Relative time within a segment;

[0387] In this embodiment, the present invention compares the sampling point of the local drone at each moment with the extrapolated position of all other drones at the same predicted moment. If the distance is lower than the preset safe distance, it is determined that there is a potential collision risk.

[0388] drones The sampling point at each moment in the predicted trajectory is compared with the position of the predicted trajectory of other drones at the same moment. If the distance is lower than the preset safe distance, it is determined that there is a potential collision risk.

[0389] When a potential collision risk is identified, each drone with this risk receives the target's status information in real time and simultaneously acquires the position information of other drones. Each drone calculates its own Euclidean distance from the target and compares it with the Euclidean distances of other drones to the target. Drones are then sorted in ascending order based on their Euclidean distances to determine their movement priority. The drone with the highest priority maintains its planned trajectory, while the remaining drones perform obstacle avoidance at different altitudes according to their priority ranking. The specific expression is:

[0390]

[0391] In the formula, For the first j The location coordinates of the drone Here are the target location coordinates, and rank(•) is the sorting function.

[0392] In the description of this invention, it should be understood that the terms "center," "thickness," "upper," "lower," "horizontal," "top," "bottom," "inner," "outer," and "radial," etc., indicating orientation or positional relationships based on the orientation or positional relationships shown in the accompanying drawings, are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of the invention. Furthermore, the terms "first," "second," and "third" are used for descriptive purposes only and should not be construed as indicating or implying the relative importance or the number of technical features implicitly specified. Therefore, a feature defined by "first," "second," and "third" may explicitly or implicitly include one or more of that feature.

Claims

1. A method for cooperative target encirclement of a UAV swarm based on brain-inspired computing, characterized in that, The method comprises the following steps: S1, acquiring multi-view images under camera view angles of the unmanned aerial vehicles, and obtaining local poses and covariances of each unmanned aerial vehicle by using panorama, visual motion and grid-landmark cell model; S2, obtaining global consistent relative poses of the unmanned aerial vehicle cluster by compression broadcasting, co-view matching and graph optimization algorithm according to the local poses and covariances of each unmanned aerial vehicle; S3, identifying the target by each unmanned aerial vehicle based on the global consistent relative poses of the unmanned aerial vehicle cluster, estimating the state of the target by using unscented Kalman filtering algorithm, and fusing all target state estimations to obtain target state information in a unified coordinate system; S4, obtaining the target motion trajectory in a historical state by using Bezier curve generation algorithm based on historical target state information, extending the Bezier curve trajectory to predict the motion of the target, and obtaining the predicted state of the target in the unified coordinate system; S5, extending the unmanned aerial vehicle cluster by using the method of centroid extension to obtain a hunting queue for hunting the target, solving the minimum hunting cost by using Gauss-Newton method to obtain an optimal hunting queue, and generating hunting positions of each unmanned aerial vehicle; S6, generating motion trajectories of the unmanned aerial vehicles by using hybrid A* search algorithm based on the hunting positions of each unmanned aerial vehicle, performing anti-collision detection according to the motion trajectories to generate a final hunting trajectory sequence, and completing hunting of the agile target; S5 comprises the following steps: S51, calculating the centroid of the unmanned aerial vehicle cluster based on the global unified coordinate; S52, calculating the vector from the centroid to each unmanned aerial vehicle, extending all the unmanned aerial vehicles outward to form a hunting queue, and sorting the unmanned aerial vehicles in the centroid area and small sector area by using permutation and combination to form several hunting queues; S53, solving the hunting function by using Gauss-Newton method to obtain the minimum hunting cost corresponding to the hunting queue; S54, comparing the minimum hunting costs corresponding to all candidate hunting queues, taking the hunting queue corresponding to the minimum minimum hunting cost as the optimal hunting queue, calculating the optimal hunting point, and obtaining the hunting position of each unmanned aerial vehicle; S52 is specifically: obtaining the vector of each unmanned aerial vehicle relative to the center of mass, calculating the angle between the vector and the surrounding direction vector of the unmanned aerial vehicle, wherein the angle between the vector and the surrounding direction vector of the first unmanned aerial vehicle is calculated as i The angle between the vector and the surrounding direction vector of the first unmanned aerial vehicle is calculated as The expression of the angle between the vector and the surrounding direction vector of the first unmanned aerial vehicle is specifically: In the formula, It is a modulo function. For vectors and x The included angle between the axes, with values ​​ranging from [0, 2π). For the first i The vector of the drone relative to its center of mass. , Y i The coordinates are in the body coordinate system. C Center of mass, , N For the number of drones, The direction vector for encirclement and capture. , T The target encirclement point; sorting all the angles of the unmanned aerial vehicles, assigning serial numbers to the unmanned aerial vehicles according to the sorting result, and generating a preliminary hunting order; if there is a unmanned aerial vehicle at the centroid position, sequentially inserting the unmanned aerial vehicle between two adjacent unmanned aerial vehicles in the circular formation, and generating a hunting queue according to all possible insertion permutation orders; if there are more than a set number of unmanned aerial vehicles in the same direction and the change of the angles of the unmanned aerial vehicles in the same sector is less than a preset threshold, arranging and combining the unmanned aerial vehicle set to obtain all possible insertion orders, inserting the unmanned aerial vehicles into the interval positions of the preliminary hunting order, and generating a hunting queue according to all possible insertion permutation orders; In S53, the minimum trapping cost corresponding to the trapping queue is calculated The expression is specifically: wherein S is a trapping queue, is an initial angular offset, R is a trapping radius, is a cost function; In the formula, to trap the queue S the middle of the j actual position of the drone.

2. The method of claim 1, wherein, S1 comprises the following steps: S11, generating a continuous panorama by using back projection and spherical mapping algorithm based on the multi-view images acquired by the unmanned aerial vehicle cluster through a monocular camera; S12, obtaining a rotation matrix and a translation vector between adjacent frames as a measurement visual motion by using ORB operator and RANSAC five-point method based on the panorama; S13, based on the feature description amplitude, using the first pulse coding and grid-landmark cell model, position coding and memory are carried out; S14, based on visual motion, IMU data and position coding, UKF is used to predict and correct the state, and the local pose and covariance of each unmanned aerial vehicle are output.

3. The method of claim 2, wherein, In S11, the method of back projection and spherical mapping algorithm is specifically: the multi-view image is back projected to the body coordinate system, wherein, for the first i Unmanned aerial vehicle camera pixels The expression of back projection is specifically: In the formula, is the coordinate of the pixel coordinate back-projection to the UAV body coordinate system, is the rotation matrix of the UAV in the world coordinate system, is the translation matrix of the starting point of the UAV in the world coordinate system, is the first i is the UAV camera intrinsic matrix, u is the horizontal direction coordinate of the image, v is the vertical direction coordinate of the image; Map all Y i On a sphere of radius to angular space and synthesize a panorama by bilinear interpolation ; S12 is specifically: according to the panoramic map, key points and their descriptors are extracted by using ORB operator, and are matched by FLANN, outliers are removed by RANSAC, and the relative motion is solved under the known intrinsic parameters by using five-point method to obtain the measured visual motion ; wherein is a rotation matrix, is a translation vector; In S13, the grid-anchor cell model comprises a grid cell and an anchor cell network, and S13 is specifically: mapping the feature description amplitude to the pulse time , to generate sparse pulses; In the formula, to prevent small amounts of division by zero, to encode the upper limit of the time window, to describe the sub, to describe the maximum amplitude at the current time; The sparse pulse is input into the grid cell and landmark cell network, and the grid cell state vector is updated by the following formula; wherein is a non-linear activation, and is a learnable parameter, is time t a grid cell state vector, is time t +1a grid cell state vector; The landmark cell activates the sparse pulse and the grid cell state vector according to the Hebbian rule to associate the memory; wherein is the activated memory of the association, is the memory before the activation of the association, is the time t is the sparse pulse, is the learning rate, is the transpose symbol; establishing correction parameters in subsequent frames based on the activated associations correcting the grid code, the corrected grid code serving as the position code wherein is time t +1 sparse pulses; In S14, the state vector of the UKF and the observation is: wherein is the time t local pose, is the time t velocity; First estimate is obtained by UKF nonlinear prediction Second estimate is obtained by EKF linear prediction Kalman gain is calculated according to second estimate ; wherein H is the state transition equation, is the noise matrix; According to the Kalman gain and the first estimation, the state vector of the UKF is updated, and the local pose and covariance of each unmanned aerial vehicle are output by the UKF; wherein is the updated state vector for the UKF, is the state transition equation transformation.

4. The method of claim 3, wherein, S2 includes the following steps: S21, based on the covariance and sparse pulse, using differential quantization or autoencoder algorithm, broadcast data packet is generated for sharing neighbor observation; S22, according to the broadcast data packet, the overlapping area of the field of view is estimated, and ORB+RANSAC algorithm is used to calculate the relative observation between the unmanned aerial vehicles; S23, based on the relative observation between the unmanned aerial vehicles, the incremental Gauss-Newton algorithm is used to solve the global pose and covariance of the unmanned aerial vehicles; S24, the global pose and covariance of the unmanned aerial vehicles are sent to each machine, and the global consistent relative pose of the unmanned aerial vehicle cluster is obtained by closed-loop correction of the global pose.

5. The method of claim 4, wherein, S21 is specifically: each unmanned aerial vehicle is different quantization of the output of the local pose increment and covariance, to get the difference quantization parameters and : wherein is a unified quantization function for vectors and matrices, is the covariance of time t is the covariance of time is the covariance of time t -1, is a local pose increment; wherein is the time t local pose of -1; The grid cell state vector and sparse pulses are reduced in dimension using a sparse autoencoder to obtain a reduced dimension result h ; wherein is a concatenation vector, is an activation function, is a sparse autoencoder parameter, is a state vector corresponding estimated value; establishing a broadcast packet , the broadcast packet arrives at the neighboring drone within a maximum latency allowed by link quality ​ S22 is specifically: the grid activation of the neighbor machine is obtained by receiving the broadcast data packet of the neighbor machine, the center deviation of the field of view is calculated according to the grid activation of the unmanned aerial vehicle itself and the grid activation of the neighbor machine, the center and radius of the overlapping circle of the field of view are estimated, and the overlapping area of the field of view is obtained; wherein the field center deviation comprises a first spherical coordinate deviation and a second spherical coordinate deviation ; wherein is the drone's own grid activation, is the neighbor drone grid activation, and is the mapping function from grid encoding to spherical coordinate difference. In the overlapping area of the field of view, the ORB operator is extracted from the decoded panoramic image segments, and the joint descriptor is constructed by combining the timing information of the first pulse coding. After RANSAC screening, the matching pairs are obtained , the relative observation between unmanned aerial vehicles is solved under the condition of known camera internal parameters K , including the rotation matrix and the translation vector of the unmanned aerial vehicle itself relative to the adjacent machine ; In the formula, is the first i frame unmanned aerial vehicle and the second j frame unmanned aerial vehicle, is the first i frame unmanned aerial vehicle and the second K frame unmanned aerial vehicle, is the first j frame unmanned aerial vehicle and the second K frame unmanned aerial vehicle; Observation covariance of relative observations between drones The expression of the observation covariance of relative observations between drones is given by: In the formula, the number of all drones that are in contact; The S23 is specifically: a graph optimization model of the UAV cluster is established, a node set is a local pose of each UAV, an edge set is composed of relative observation and observation covariance between the UAVs, and an optimization objective function of the graph optimization model is The expression of the optimization objective function is specifically: wherein is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, j is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, i is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, is the pose change difference between the initial position and the position of the third UAV in the three-dimensional space, is the SE(3) composition; The working process of the incremental Gauss-Newton algorithm is as follows: when receiving new relative observations, the residual error is linearized , the sparse information matrix is incrementally updated, the local pose of the UAV after updating is solved, and the global pose of the UAV is obtained. In the formula, For the objective function The first partial derivative, For the objective function The first differential, For the objective function The second partial derivative, For the objective function The second differential, This is the difference between the sum of the partial derivative and the differential product; In the formula, is the sparse information matrix before update, is the sparse information matrix after update; In the formula, p i to update the local pose of the UAV before, to update the local pose of the UAV after S24 is specifically: the global pose and covariance of the unmanned aerial vehicle are broadcast to each unmanned aerial vehicle, and the UKF state update of the local unmanned aerial vehicle is specifically: wherein, is the updated state vector, is the global pose of the UAV, is the updated grid cell state vector, P i is the covariance, and is the prior noise for velocity and grid activation, is the updated pose in three-dimensional space of the UAV, is the updated velocity; Each drone utilizes and restarts the filtering loop and updates the long-term memory of the updated grid cell state vector and sparse pulse input landmark cells to obtain the global consistent relative pose of the drone swarm.

6. The method of claim 5, wherein, S3 includes the following steps: S31, using YOLO target recognition method, the target in the image information is recognized, and the target position in the coordinate system thereof is obtained according to the image information; S32, the target position is transformed through the rotation matrix to obtain the target position in the unified coordinate system; S33, based on the target position in the unified coordinate system, UKF is used to estimate the state of the target, and the target state estimation and covariance after filtering the sensor error are obtained; S34, in a loosely coupled manner, all target state estimations are fused to obtain target state information by combining the observation noise covariance of each unmanned aerial vehicle on the target; S33 is specifically: The UKF generates 2 +1 sigma points from the state vector and covariance matrix n according to the following formula: wherein is a sigma point, n is the dimension of the state vector, is K is the estimate of the position at the time instant, is a tuning function; wherein and are tuning parameters; propagating the sigma points through a nonlinear state transition function to compute a predicted state mean and covariance ; wherein is obtained by K the state value prediction at time K the value at time Q is a covariance matrix of the measurement noise, , is the state part, is the noise part, and is a weight coefficient; In the formula, is a small error; mapping the predicted Sigma points to the measurement space, calculating the mean of the measurements , the covariance matrix S and the cross-covariance matrix T ; wherein is obtained by K the observation at time K the value at time R i is a predetermined noise covariance matrix assumed Calculating Kalman gain from covariance matrix ; According to the observation value and the measured mean value The state estimation and covariance of the target are updated, and the updated state estimation of the target and covariance The expression is: In the formula, is K an observation of the time In S34, the target state information is obtained The expression is specifically: In the formula, is the first i target state estimation of the UAV, is the first weight of the UAV; In the formula, is the first The observation noise covariance of the UAV is obtained.

7. The method of claim 6, wherein, S4 includes the following steps: S41, using sliding window method to record historical target state information with a set length; S42, based on the historical target state information, a target motion trajectory in the historical state is generated according to the Bezier curve generation algorithm; S43, according to the distance between the centroid of the unmanned aerial vehicle cluster and the current target position, the predicted step is calculated by hyperbolic tangent function; S44, the Bezier curve is extended by the predicted step to obtain the predicted state of the target in the unified coordinate; S42 is specifically: using a Bezier curve to describe the target motion trajectory, and the expression of the target motion trajectory is specifically: wherein is an n-th Bernstein polynomial basis, are Bezier curve control points; In S43, the predicted step size is calculated The expression of the step size is specifically: In the formula, is a preset prediction time step, d is a target current position distance, is a current time.

8. The method of claim 7, wherein, S6 includes the following steps: S61, obtain surrounding obstacle information from binocular camera, construct three-dimensional grid or octree environment model, map continuous state of UAV to discrete grid and heading discrete angle, for search on grid map, continuous state includes position, orientation and velocity; S62, take current position, orientation and velocity of UAV as search starting point, determine safe surrounding point and posture of target region or target swarm according to predetermined "capture enclosure" strategy; S63, in discretized state space, take kinematic constraint of UAV as expansion rule, use heuristic function to evaluate cost of current grid point to target, push search to develop towards target direction, for generating trajectory satisfying dynamics constraint and being connected; S64, perform spline fitting or smoothing processing of minimum bending energy on discrete trajectory generated based on hybrid A* search algorithm, eliminate sharp turns of polyline segment, assign velocity and time stamp to smooth trajectory according to velocity and acceleration limit of UAV, obtain continuous executable time parameterized trajectory; S65, discretely sample smooth trajectory, use nearest obstacle distance detection algorithm to detect obstacle collision at each sampling point, in response to detecting collision or insufficient safety distance, then trigger local re-planning, perform high-low flight in the vicinity of the segment to correct path in smooth trajectory; S66, integrate and splice smooth trajectory, time parameterized trajectory and trajectory segment that has passed collision detection, generate final capture trajectory sequence; S64 is specifically: Target-oriented kinematic search front-end: search tree is constructed based on hybrid A* search algorithm, motion primitives are generated by discretizing control inputs, and cost function is introduced in search process With target trajectory B(t) as heuristic information, double heuristic function is designed: wherein, is a dynamic distance estimation based on optimal boundary value problem, is a UAV target state, is a UAV current state, is a weight, S t is a sum of expected expansion times, is a path expansion time, is obtained by extrapolation of the predicted trajectory: In the formula, a Bezier curve trajectory of the original historical trajectory, a trajectory obtained by time extrapolation; spatial trajectory optimization backend: constructs a flight corridor consisting of connected free-space cubes based on an initial path generated by a goal-oriented kinematic search frontend ; wherein is a Bezier curve control point, M is the number of trajectory segments of the flight corridor; Establish optimization problem with dynamics constraint under corridor constraint: wherein is a smoothing term, is a potential field term, is a dynamic constraint term, q w is an intermediate waypoint, T is a time of the piecewise trajectory, is x , y and z a combination of, x is an x axis coordinate in the piecewise waypoint, y is an y axis coordinate in the piecewise waypoint, z is an z axis coordinate in the piecewise waypoint; the optimization problem is solved efficiently by a split convex optimization method, generating a time-parameterized trajectory; In S66, the method of collision detection is specifically: For each time-parameterized trajectory of a drone, let the trajectory of the drone be composed of piecewise polynomials, define a set of parameters m ;​ wherein is the k segment duration, , is the real domain, is the k segment polynomial order, , is the matrix domain, is the coefficient matrix, , is the k segment polynomial's x parameter in the axis direction, is the k segment polynomial's y parameter in the axis direction, is the k segment polynomial's z parameter in the axis direction; Calculating the relative time in each segment at the current time according to the polynomial ; wherein, is the difference between the current time and the start time, is the absolute time, is the track start time stamp, is the max function, s is the number of track segments, is the min index function, is the time of each track segment; Computing a three-dimensional spatial position of a predicted trajectory ; In the formula, is the current position of the UAV, is the first s segment, j is the coefficient of the , is the x coefficient in the axis direction, is the y coefficient in the axis direction, is the z coefficient in the axis direction, is the relative time in the first s segment. The unmanned aerial vehicle The sampling points of each time in the predicted trajectory are compared with the positions of other unmanned aerial vehicles at the same time in the predicted trajectories, and if the distance is lower than a preset safety distance, it is determined that there is a potential collision risk. When it is determined that there is a potential collision risk, each unmanned aerial vehicle in the potential collision risk receives state information of the target in real time, and acquires position information of other unmanned aerial vehicles, each unmanned aerial vehicle calculates the Euclidean distance between itself and the target, compares the Euclidean distance with the Euclidean distance between other unmanned aerial vehicles and the target, sorts the Euclidean distances in ascending order according to the size of the Euclidean distance to determine the motion priority, the unmanned aerial vehicle with the highest priority keeps moving according to the original planned trajectory, and the remaining unmanned aerial vehicles avoid obstacles at different heights according to the priority sorting, and the unmanned aerial vehicles avoid obstacles at different heights according to the sequence number after the priority sorting The expression is specifically: In the formula, is the first j The position coordinates of the UAV on the shelf, is the target position coordinates, rank( ) is the ranking function

Images