Non-cooperative UAV Trajectory Distribution Prediction Method Based on Flight State Division

Abstract

The present invention discloses a non-cooperative UAV trajectory distribution prediction method based on flight state division, and the steps are as follows: for the selected monitoring airspace range, construct a rasterized airspace; divide the UAV flight states; screen the similar trajectory data sets and construct a trajectory prediction model based on data migration; generate the trajectory reachable space; index the grid coordinates covered by the trajectory; generate the trajectory probability distribution. The method of the present invention realizes more accurate trajectory distribution prediction and reduces the range of the trajectory space with non-zero probability by subdividing the flight states of the UAV, combining the trajectory prediction method based on data migration, modeling the motion considering the uncertainty of the UAV operator's intention as Brownian motion, and using the truncated Brownian bridge method to model the position distribution of non-cooperative UAVs.

Description

Technical Field

[0001] The present invention belongs to the technical field of unmanned aerial vehicle (UAV) traffic management, and particularly relates to a method for predicting the trajectory distribution of non-cooperative UAVs based on flight state classification. Background Art

[0002] With the rapid popularization of UAVs in various industries, the development of low-altitude UAVs is rapid and the flight demand is increasing day by day, and air traffic activities are becoming increasingly frequent. The research on UAV trajectory prediction technology can promote the safe, orderly, and efficient operation of UAVs in the limited low-altitude airspace. Calculating the trajectory distribution of UAVs at a certain moment in the airspace is the basis for constructing a dangerous behavior recognition model.

[0003] Most of the existing trajectory prediction methods obtain the estimated position points of UAVs in the future, rather than the position intervals where the UAVs are located. When analyzing the safety risks of UAV operations (judging the intrusion of a fixed area or conflict risk), it is possible to ignore the original dangerous behaviors or generate false alarms due to the errors in the position prediction results. In addition, most of the existing trajectory distribution predictions of UAVs are based on the maximum reachable range of the trajectory or the reachable range of the trajectory under some known conditions. However, in the actual task scenarios that require real-time trajectory prediction, it may face non-cooperative targets that appear for the first time, and their identity information is often unknown during the implementation of monitoring and management. It is very difficult to model and estimate the intentions of the UAV operators of this type. Summary of the Invention

[0004] Aiming at the deficiencies of the above-mentioned existing technologies, the purpose of the present invention is to provide a method for predicting the trajectory distribution of non-cooperative UAVs based on flight state classification, so as to solve the problems that it is difficult to estimate the intentions of UAV operators when the existing trajectory prediction methods face non-cooperative targets, and the prediction accuracy of future trajectory points is not high, resulting in missed alarms or false alarms. The method of the present invention divides the initial flight state of the UAV, combines probability distribution prediction to replace the prediction of future trajectory points, generates a trajectory distribution space, reduces the occurrence of missed alarms and false alarms, and at the same time reduces the range of the trajectory space with non-zero probability, realizing a more accurate trajectory distribution prediction.

[0005] To achieve the above purpose, the technical solution adopted by the present invention is as follows:

[0006] A method for predicting the trajectory distribution of non-cooperative UAVs based on flight state classification of the present invention comprises the following steps:

[0007] (1) For the range of the selected monitoring airspace, construct a rasterized airspace;

[0008] (2) Divide the flight state of the UAV: According to the flight state of the UAV at the initial moment, divide the flight state of the UAV into a hovering state and a moving state, and set the UAV trajectory distribution prediction parameters;

[0009] (3) Screen the similar trajectory datasets and construct a trajectory prediction model based on data migration;

[0010] (4) Generate the trajectory reachable space: According to the flight states of the UAVs divided in step (2), solve the set of all possible trajectory positions of the UAVs at the prediction moment as a region with a closed boundary, and generate the trajectory reachable space of the UAVs at the prediction moment;

[0011] (5) Index the grid coordinates covered by the trajectory: Approximately express the trajectory reachable space generated in step (4) as the set of grids it covers in the rasterized airspace, and solve the coordinates of the covered grids in the Cartesian coordinate system of the airspace;

[0012] (6) Generate the trajectory probability distribution: Solve the probability that the UAV is located in each grid in the Cartesian coordinate system of the airspace at the monitoring moment, and obtain the probability distribution corresponding to the set of grids.

[0013] Furthermore, for the monitoring airspace range in step (1), establish a Cartesian coordinate system LL of the airspace with the lower left corner as the origin. By setting the number of grids on the x-axis and y-axis, divide the airspace into a set of grids with unique coordinates. The grid coordinates (i, j) indicate that the grid is the i-th grid in the x-axis direction and the j-th grid in the y-axis direction of the airspace. The information saved in each grid includes: the coordinates of the grid center point and the four vertexes in the Cartesian coordinate system LL of the airspace.

[0014] Furthermore, the specific process of step (2) is as follows:

[0015] (21) Divide the flight states of the UAVs;

[0016] According to the flight states of the UAVs at the initial moment, divide the flight states of the UAVs into a hovering state and a moving state;

[0017] (22) Set the prediction parameters of the UAV trajectory distribution;

[0018] (221) Set the initial flight speed of the UAV;

[0019] For the initial speed of the UAV during the trajectory prediction process, set the absolute value of the initial speed in the hovering state to zero, and the absolute value of the initial speed in the moving state to be greater than zero;

[0020] (222) Set the maximum ground speed and maximum horizontal acceleration of the UAV;

[0021] According to the factory performance parameters of the UAV, estimate that the maximum allowable flight speed of the UAV under calm wind is less than 30 m / s, set the maximum ground speed to 30 m / s, and the maximum horizontal acceleration to 6 m / s 2 ;

[0022] (223) Set the prediction duration of the UAV trajectory;

[0023] According to the maximum flight ground speed, maximum horizontal acceleration, maximum static wind speed of the UAV set in step (222) and the pre-time requirement for UAV dangerous behavior recognition, set the prediction duration to be greater than the time required for the UAV to accelerate from the hovering state to the maximum flight ground speed at the maximum horizontal acceleration.

[0024] Furthermore, the specific process of step (3) is as follows:

[0025] (31) Screen the similar trajectory dataset;

[0026] Collect the trajectory segments of non-cooperative UAVs. Based on the existing cooperative UAV trajectory dataset, express the potential characteristics of the UAV trajectory through the flight speed variance and the cumulative heading change amount, and then describe the stability of the flight state of each trajectory segment. Measure the distance between trajectories based on the similarity of the flight speed variance and the cumulative heading change amount values, and screen out the trajectory dataset similar to the non-cooperative UAV trajectory;

[0027] (32) Construct a trajectory prediction model based on data migration;

[0028] Use the similar trajectory dataset screened in step (31) as the training sample of the D-GRU trajectory prediction model, and regard the trained D-GRU trajectory prediction model as the trajectory prediction model based on data migration.

[0029] Furthermore, the specific process of step (4) is as follows:

[0030] (41) Spatiotemporal constrained Brownian bridge;

[0031] For the motion characteristics of the UAV in different flight states divided in step (21), when representing the position of the UAV at time t as (x(t), y(t)) according to the truncated Brownian bridge, x(t) follows the truncated probability density function:

[0032]

[0033] At time t, the UAV reaches the coordinate in the trajectory reachable space For any y(t) follows the truncated probability density function:

[0034]

[0035] In addition, the probability that the UAV reaches the coordinate in the trajectory reachable space is:

[0036]

[0037] In the formula, (x(t), y(t)) is the position of the UAV at time t. For the normal distribution X~N(μ, σ 2 ), is the probability density function, Φ(x) is the distribution function, U x (t), U y (t) are the upper boundaries in the x and y directions, and L x (t), L y (t) are the lower boundaries in the x and y directions;

[0038] (42) Reachable space of the UAV trajectory in the initial hovering state;

[0039] For a UAV in the hovering state at the initial moment, the reachable space of the trajectory at the prediction moment is a circle with the center at the position at the initial moment and the radius being the maximum flight distance to the prediction moment. The maximum flight distance includes two parts: the accelerated flight distance and the uniform flight distance at the maximum flight speed;

[0040] (43) Reachable space of the UAV trajectory in the initial motion state;

[0041] (431) Determine the minimum turning radius;

[0042] The minimum turning radius of the UAV is obtained based on the maximum overload factor allowed during the design of the UAV:

[0043]

[0044] In the formula, R min is the minimum turning radius of the UAV, g represents the acceleration due to gravity, |v 0 | represents the magnitude of the flight speed of the UAV when starting to turn, and n max is the maximum overload factor. It is assumed that the UAV with an initial speed starts to turn at the initial moment, and the initial speed remains unchanged in magnitude during the turn, only changing the direction;

[0045] (432) Set the direction of the turning angle;

[0046] The difference between the UAV speed directions before and after turning is defined as the turning angle. When turning to the right, the turning angle is positive; when turning to the left, the turning angle is negative; when the initial speed direction remains unchanged, the turning angle is 0;

[0047] (433) Establish the Cartesian coordinate system of the UAV;

[0048] Taking the position of the UAV at the initial moment as the coordinate origin, the direction of the initial speed of the UAV as the positive direction of the vertical axis Y, and the right side of the fuselage as the positive direction of the horizontal axis X, establish the Cartesian coordinate system OXY of the UAV;

[0049] (434) Analyze the reachable space of the trajectory according to different motion phases;

[0050] Based on the coordinate system established in step (433), for the turning phase of the UAV, solve the coordinates of the UAV after the turn according to the direction of the turning angle in step (432); for the acceleration phase of the UAV, set the acceleration from the initial speed to the maximum flight speed with the maximum horizontal acceleration; for the uniform flight phase of the UAV, fly at a constant speed with the maximum flight speed; solve the position of the UAV at the prediction moment, that is, the last position in the uniform flight phase;

[0051] (435) Discretize the boundary of the reachable space of the trajectory;

[0052] The boundary of the reachable space of the UAV's trajectory at the prediction moment is generated by the set of boundary coordinates of the UAV at that moment with the turning angle in the interval (-360, 360). Discretize the turning angle, sample the turning angle at equal interval degrees, and use the end point of the previous flight phase as the starting point of the next phase to generate the reachable space of the trajectory;

[0053] (436) Clear the overlapping areas of the reachable space of the trajectory;

[0054] Solve the magnitude of the turning angle of the acceleration trajectory when the first intersection appears in the flight trajectory of the UAV during the time lengths of turning and accelerating operations:

[0055]

[0056] In the formula, s a is the flight distance of the UAV during the acceleration phase, and R min is the minimum turning radius of the UAV; clear the reachable space of the trajectory where overlap occurs during large-angle turns;

[0057] (437) Regenerate the reachable space of the trajectory for different motion phases;

[0058] Sample the turning angle interval obtained in step (436) at equal interval degrees to generate the reachable space of the trajectory covered in different phases of the initial motion.

[0059] Furthermore, the specific process of step (5) is as follows:

[0060] (51) Index the grid coordinates covered by the trajectory in the initial hovering state;

[0061] (511) Generate the initial coverage space;

[0062] Generate the initial coverage space according to the circumscribed square of the circle in step (42), and save the coordinate sets of the upper left, lower left, upper right, and lower right vertices of the circumscribed square;

[0063] (512) Expand the initial coverage space;

[0064] Based on the possible incomplete grid coverage at the boundary of the initial coverage space in step (511), for the grid coordinates of the lower boundaries in the x and y directions, round down, and for the grid coordinates of the upper boundaries, round up to obtain the expanded coverage space;

[0065] (513) Solve for the coordinate set of the grid set in the gridified airspace Cartesian coordinate system LL;

[0066] Traverse the grids included in the coverage space in step (512), and retain the grids whose distance from the grid center to the initial position of the UAV is less than or equal to the maximum flight distance to obtain the coordinate set of the grid set in the gridified airspace Cartesian coordinate system LL;

[0067] (52) Initial motion state trajectory coverage grid coordinate index;

[0068] (521) Construct a discrete boundary coordinate set;

[0069] Based on the motion state of the UAV at the initial moment, obtain the discrete boundary coordinate set of the reachable space of its trajectory at the predicted moment in the UAV Cartesian coordinate system OXY in step (433);

[0070] (522) Coordinate system transformation;

[0071] For each discrete boundary coordinate (x c，β , y c，β ) in the UAV Cartesian coordinate system OXY in step (521), convert it to the coordinate (x L c，β , y L c，β ) in the airspace Cartesian coordinate system LL where the grid is located in step (1):

[0072]

[0073] In the formula, α L，o represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. The initial position coordinates of the UAV are (x o , y o ) to obtain the boundary coordinate set on the airspace Cartesian coordinate system LL;

[0074] (523) Obtain the coverage space;

[0075] Generate the initial coverage space according to step (511), and then expand the initial coverage space according to step (512), that is, obtain the coverage space included in the minimum integerized grid of the reachable space of the trajectory;

[0076] (524) Curve fitting circle equation;

[0077] Based on the trajectory reachable space obtained in step (43), transfer the grid center coordinates (x L ck , y L ck ) to (x ck , y ck ) in the UAV Cartesian coordinate system OXY:

[0078]

[0079] In the formula, α L，O represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. The initial position coordinates of the UAV are (x o , y o ). Furthermore, calculate the angle β ck between the line connecting the point (x ck , y r，o ) in the UAV Cartesian coordinate system OXY and the initial position O(0, 0) of the UAV and the positive Y-axis:

[0080]

[0081] In the formula, (x ck , y ck ) are the coordinates of the center of any grid ck in the UAV Cartesian coordinate system OXY, β r，o ∈[-180, 180]. Furthermore, construct an arc fitting model fit by the least squares method:

[0082] (x ck，c , y ck，c ) = fit({x c，β , y c，β}), β ∈ [β r，o - ε, β r，o + ε] (9)

[0083] The input of the model is the set of boundary coordinates X k，r = {(x c，β , y c，β )}, β ∈ [β r，o - ε, β r，o + ε], where ε is used to define the angle range. Consider the points within the small angle range in the set of boundary coordinates X k，r = {(x c，β , y c，β )} as an arc, and perform curve fitting by the least squares method. The output is the center (x ck，c , y ck，c ) of the arc and the radius r ck，c, obtain the fitting circle equation corresponding to the minimum error;

[0084] (525) Determine the positional relationship between the grid center coordinates and the fitting circle;

[0085] Calculate the Euclidean distance between the grid center coordinates and the center of the fitting circle in step (524). By comparing the magnitude of this Euclidean distance with the radius of the fitting circle, determine whether the grid center coordinates are located inside the fitting circular arc; regard the grids located inside the circular arc as being within the reachable space of the trajectory;

[0086] (526) Solve the coordinate set of the grid set in the grid - based airspace Cartesian coordinate system LL;

[0087] Traverse the grids contained in the coverage space in step (523), and determine the positional relationship between the grid center coordinates and the fitting circle through step (525) to obtain the coordinate set of the grid set in the grid - based airspace Cartesian coordinate system LL;

[0088] (53) Represent the indexed grid coordinates;

[0089] Represent the coordinate sets obtained by indexing the grid sets in different flight states as the set {ck}, where ck=(i k , j k ) indicates that the grid ck in the grid set is the i - th k , j - th k grid in the x - direction and y - direction of the airspace Cartesian coordinate system LL.

[0090] Further, the specific process of step (6) is as follows:

[0091] (61) Solve the trajectory probability distribution for the initial hovering state;

[0092] (611) Coordinate system translation;

[0093] Consider the subsequent movement of the unmanned aerial vehicle (UAV) in the hovering state at the initial moment as an undirected Brownian motion, and translate the reachable space of the trajectory from the original (x, y) in the airspace Cartesian coordinate system LL to a new coordinate system (x′, y′) with the position of the UAV at the initial moment as the origin:

[0094]

[0095] where (x 0 , y 0 ) is the position coordinate of the UAV at the initial moment;

[0096] (612) Solve the trajectory probability distribution;

[0097] In the coordinate system translated in step (611), obtain the predicted - moment UAV position coordinates through formula (3) The probability density function, convert the grid coordinates in the grid set of step (513) through formula (10), grid ck = (i k , j k ). The coordinates of the grid center, lower left, upper left, upper right, and lower right vertices of the grid in the translation coordinate system are successively {(x ck ′(t), y ck ′(t)), (x ck，l ′(t), y ck，l ′(t)), (x ck，l ′(t), y ck，u ′(t)), (x ck，r ′(t), y ck，u ′(t)), (x ck，r ′(t), y ck，l ′(t))}, solve the position of the UAV in the reachable space of the trajectory at the prediction moment, i.e., at time t The probability p 1，ck (t) of being located within grid ck:

[0098]

[0099] In the formula, x ck，l ′(t), x ck，r ′(t) are the abscissas of the upper left and lower left, upper right and lower right vertices of the grid respectively, and y ck，l ′(t), y ck，u ′(t) are the ordinates of the lower left and lower right, upper left and upper right vertices of the grid respectively, is the probability density function of the coordinates at time t in the initial hovering state ;

[0100] (613) Update the trajectory probability distribution;

[0101] Sum the probability values in each grid obtained according to step (612) to obtain the total probability p 1 (t), p 1 (t) is extremely close to 1 when the grid is fine enough:

[0102]

[0103] In the formula, p 1，ck (t) is the probability that the position of the UAV in the reachable space of the trajectory at time t in the initial hovering state is located within grid ck. To extend the model to scenarios with a larger grid granularity, update p 1，ck (t) through formula (12):

[0104]

[0105] When obtaining the initial hovering state of the UAV, the probability set of the UAV at each grid in the grid set at time t is obtained;

[0106] (62) Solve the probability distribution of the initial motion state trajectory;

[0107] (621) Solve the expectation of the reachable space of the trajectory at the prediction time;

[0108] The subsequent motion of the UAV with the initial motion state is regarded as the combination of undirected Brownian motion and directed Brownian motion. Use the trajectory prediction model based on data migration in step (32) to solve the trajectory prediction value of the UAV at the prediction time, and regard it as the expectation E=(μ x (t), μ y (t));

[0109] (622) Translate and rotate the coordinate system;

[0110] Translate the reachable space of the trajectory from the grid-based airspace Cartesian coordinate system LL to the origin with the expected value E of the UAV position at the prediction time obtained in step (621), and rotate it to the positive X-axis of the new coordinate system EXY with the line connecting the coordinates at the initial time and E. Perform coordinate transformation of the original coordinate system in the new coordinate system:

[0111]

[0112] In the formula, α L，E represents the angle between the airspace Cartesian coordinate system LL and the positive X-axis of the EXY coordinate system;

[0113] (623) Determine the upper and lower limits of the abscissa of the UAV in the reachable space of the trajectory at the prediction time;

[0114] Convert the boundary coordinate set on the airspace Cartesian coordinate system LL to the coordinate system EXY in step (522), and represent the upper and lower limits of the abscissa of the UAV in the reachable space of the trajectory at the prediction time;

[0115] (624) Determine the upper and lower limits of the ordinate based on the abscissa;

[0116] Calculate the absolute value of the difference between the abscissa in step (623) and the converted boundary coordinate set, and solve the index i corresponding to the minimum value in this set x，max and the y corresponding to this index i , if y i ≥0, set the upper limit of the ordinate as y i , execute the loop, set the set index i x，max at this position as +∞, and solve the index i corresponding to the minimum value in this set x， m ax , until the index i x，max corresponding yi < 0, stop the loop and set the lower limit of the vertical coordinate to y i ; If y i < 0, set the lower limit of the vertical coordinate to y i , execute the loop, set the index i of this set x，max to +∞, and solve for the index i corresponding to the minimum value of this set x，max , until the index i x，max corresponding y i ≥0, stop the loop and set the upper limit of the vertical coordinate to y i ;

[0117] (625) Solve the probability density function at the prediction time;

[0118] According to the upper and lower limits of the abscissa determined in step (623) and the upper and lower limits of the vertical coordinate determined in step (624), obtain the probability density function of the UAV position coordinates at the prediction time when the initial state is in motion through formula (3) within the reachable space of the trajectory;

[0119] (626) Select the integration range of the probability density function;

[0120] Based on the vertex coordinates of the rotated grid in step (622), the binary integration domain will become larger. Calculate the area s of the true grid integration domain obtained based on the rotation angle of the coordinate system 1 , establish the true grid integration domain s ck and the area s of the binary integration domain of the rotated grid 2 The relationship between them is:

[0121]

[0122] (627) Solve the trajectory probability distribution;

[0123] In the translated and rotated coordinate system of step (622), convert the grid coordinates in the grid set of step (526) through formula (14). The grid center, lower left, upper left, upper right, and lower right vertices of the grid ck=(i k , j k ) in the translated coordinate system are {(x ck ′(t), y ck ′(t)), (x ck，ll ′(t), y ck，ll ′(t)), (x ck，lu ′(t), y ck，lu ′(t)), (x ck，ru ′(t), y ck，ru ′(t)), (x ck，rl ′(t), y ck，rl′(t))}, and based on the relationship obtained in step (626), establish the double integrals p 1 and p 2 on areas s 2，ck1 (t) and p 2，ck2 (t) respectively:

[0124]

[0125]

[0126] where (x ck ′(t), y ck ′(t)) is the central coordinate of grid ck, a and b are respectively half of the lengths of two adjacent sides in area s 1 , U ck，x′ (t) and L ck，x′ (t) are respectively the maximum and minimum values of the abscissas of the four vertices (lower left, upper left, upper right, lower right) of grid ck, and U ck，y′ (t) and L ck，y′ (t) are respectively the maximum and minimum values of the ordinates of the four vertices (lower left, upper left, upper right, lower right) of grid ck. is the probability density function of the coordinate at the initial motion state at time t, and according to formula (15), the position of the UAV within the reachable space of the trajectory at the initial motion state at time t is located within grid ck with probability p (t): 2，ck (t):

[0127]

[0128] (628) Solve the probability distribution set;

[0129] Update the probability p 2，ck (t) in step (627) according to the probability update formula (13) in step (613) to obtain the probability set of the UAV located in each grid in the grid set at the initial motion state of the UAV.

[0130] Advantages of the present invention:

[0131] The method of the present invention improves the accuracy of predicting the trajectory distribution of the UAV when the trajectory information is unknown, provides a basis for identifying dangerous behaviors and assessing conflict risks during the flight of the UAV, and provides theoretical support for low-altitude air traffic safety management.

[0132] The method of the present invention realizes more accurate trajectory distribution prediction and reduces the trajectory space range of non-zero probability by subdividing the flight state of the UAV, combining the trajectory prediction method based on data migration, modeling the motion considering the uncertainty of the UAV operator's intention as Brownian motion, and using the truncated Brownian bridge method to model the position distribution of non-cooperative UAVs. BRIEF DESCRIPTION OF THE DRAWINGS

[0133] Figure 1 It is a flowchart of the method of the present invention.

[0134] Figure 2a It is a schematic diagram of the trajectory reachable space of the initial hovering state in the embodiment when the prediction duration is 10 s.

[0135] Figure 2b It is a schematic diagram of the trajectory reachable space of the initial hovering state in the embodiment when the prediction duration is 20 s.

[0136] Figure 2c It is a schematic diagram of the trajectory reachable space of the initial hovering state in the embodiment when the prediction duration is 30 s.

[0137] Figure 3a It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the turning stage when the prediction duration is 5 s.

[0138] Figure 3b It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the acceleration stage when the prediction duration is 5 s.

[0139] Figure 3c It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the uniform motion stage when the prediction duration is 5 s.

[0140] Figure 4a It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the turning stage when the prediction duration is 30 s.

[0141] Figure 4b It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the acceleration stage when the prediction duration is 30 s.

[0142] Figure 4c It is a schematic diagram of the trajectory reachable space of the initial motion state in the embodiment at the uniform motion stage when the prediction duration is 30 s.

[0143] Figure 5 It is a schematic diagram of the trajectory reachable space S p,t , the initial coverage space S p,t,1 and the minimum integerized grid S p,t,2 relationship diagram.

[0144] Figure 6a Schematic diagram of grid integration interval when the coordinate system rotation angle is 12° in the embodiment of the present invention.

[0145] Figure 6b Schematic diagram of grid integration interval when the coordinate system rotation angle is 50° in the embodiment of the present invention. Detailed implementation manners

[0146] For the convenience of those skilled in the art to understand, the present invention will be further described below in conjunction with embodiments and the accompanying drawings. The content mentioned in the implementation manners does not limit the present invention.

[0147] Refer to Figure 1 As shown, a non-cooperative UAV trajectory distribution prediction method based on flight state division of the present invention is as follows:

[0148] (1) For the selected monitoring airspace range, construct a grid-based airspace;

[0149] Among them, the monitoring airspace range establishes an airspace Cartesian coordinate system LL with the lower left corner as the origin. By setting the number of grid cells on the x-axis and y-axis, the airspace is divided into a grid set with unique coordinates. The grid coordinate (i, j) indicates that the grid is the i-th grid in the x-axis direction and the j-th grid in the y-axis direction of the airspace. The information stored in each grid includes: the coordinates of the grid center point and the four vertex points in the airspace Cartesian coordinate system LL.

[0150] (2) Divide the UAV flight states: According to the flight state of the UAV at the initial moment, divide the flight state of the UAV into a hovering state and a moving state, and set the UAV trajectory distribution prediction parameters;

[0151] (21) Divide the UAV flight states;

[0152] According to the flight state of the UAV at the initial moment, divide the flight state of the UAV into a hovering state and a moving state;

[0153] (22) Set the UAV trajectory distribution prediction parameters;

[0154] (221) Set the initial UAV flight speed;

[0155] For the initial speed of the UAV during the trajectory prediction process, set the absolute value of the initial speed in the hovering state to zero, and the absolute value of the initial speed in the moving state to be greater than zero;

[0156] (222) Set the maximum UAV flight ground speed and the maximum horizontal acceleration;

[0157] According to the factory performance parameters of the UAV, it is estimated that the maximum allowable flight speed of the UAV under calm wind is less than 30 m / s. Set the maximum ground speed to 30 m / s and the maximum horizontal acceleration to 6 m / s 2 ;

[0158] (223) Set the UAV trajectory prediction duration;

[0159] According to the maximum ground speed, maximum horizontal acceleration, maximum calm wind speed of the UAV set in step (222) and the pre-time requirement for UAV dangerous behavior recognition, set the prediction duration to be greater than the time required for the UAV to accelerate from the hover state to the maximum ground speed at the maximum horizontal acceleration.

[0160] (3) Screen the similar trajectory dataset and construct a trajectory prediction model based on data migration;

[0161] (31) Screen the similar trajectory dataset;

[0162] Collect the trajectory segments of non-cooperative UAVs. According to the existing cooperative UAV trajectory dataset, express the potential characteristics of the UAV trajectory through the flight speed variance and the cumulative heading change amount, and then describe the stability of the flight state of each trajectory segment. Based on the similarity of the flight speed variance and the cumulative heading change amount values, measure the distance between trajectories and screen out the trajectory dataset similar to the non-cooperative UAV trajectory;

[0163] (32) Construct a trajectory prediction model based on data migration;

[0164] Use the similar trajectory dataset screened in step (31) as the training sample of the D-GRU trajectory prediction model, and regard the trained D-GRU trajectory prediction model as the trajectory prediction model based on data migration.

[0165] (4) Generate the trajectory reachable space: According to the UAV flight states divided in step (2), solve the set of all possible trajectory positions of the UAV at the prediction moment as a region with a closed boundary, and generate the trajectory reachable space of the UAV at the prediction moment;

[0166] (41) Spatiotemporal constrained Brownian bridge;

[0167] For the motion characteristics of the UAV in different flight states divided in step (21), when expressing the position of the UAV at time t as (x(t), y(t)) according to the truncated Brownian bridge, x(t) follows the truncated probability density function:

[0168]

[0169] The UAV reaches the coordinate within the trajectory reachable space at time t For any y(t) follows a truncated probability density function:

[0170]

[0171] In addition, the probability that the UAV reaches the coordinate within the reachable space of the trajectory is:

[0172]

[0173] where (x(t), y(t)) is the position of the UAV at time t. For a normal distribution X ∼ N(μ, σ 2 ), is the probability density function, Φ(x) is the distribution function, U x (t), U y (t) are the upper boundaries in the x and y directions, and L x (t), L y (t) are the lower boundaries in the x and y directions;

[0174] (42) Reachable space of the UAV trajectory in the initial hovering state;

[0175] For a UAV initially in a hovering state, the reachable space of the trajectory at the prediction time is a circle with the center at the position at the initial time and the radius equal to the maximum flight distance to the prediction time. The maximum flight distance includes two parts: the accelerating flight distance and the constant-speed flight distance at the maximum flight speed. The prediction durations are taken as 10 s, 20 s, and 30 s. The reachable space of the UAV trajectory sampled at 60-degree intervals is as shown in Figures 2a - 2c . The lines shown in the figure are the headings of 6 UAVs sampled at 60-degree intervals. The thinner lines closer to the center are the accelerating flight stages of the UAVs, and the thicker lines closer to the outside are the stages of the UAVs flying at a constant speed at the maximum flight speed. As the prediction duration of dangerous behavior identification increases, the proportion of the constant-speed flight distance at the maximum flight speed in the radius will increase significantly;

[0176] (43) Reachable space of the UAV trajectory in the initial motion state;

[0177] (431) Determine the minimum turning radius;

[0178] The minimum turning radius of the UAV is obtained based on the maximum overload factor allowed during the design of the UAV:

[0179]

[0180] where R min is the minimum turning radius of the UAV, g represents the acceleration due to gravity, |v 0 | = 10 m / s represents the magnitude of the flight speed of the UAV when starting to turn, and n max= 3.5 is the maximum overload factor. It is assumed that a UAV with an initial velocity starts to turn at the initial moment. During the turn, the magnitude of the initial velocity remains unchanged while only the direction changes;

[0181] (432) Set the direction of the turning angle;

[0182] Define the difference between the velocity directions of the UAV before and after turning as the turning angle. When turning right, the turning angle is positive; when turning left, the turning angle is negative; when the initial velocity direction remains unchanged, the turning angle is 0;

[0183] (433) Establish a Cartesian coordinate system for the UAV;

[0184] Taking the position of the UAV at the initial moment as the coordinate origin, the direction of the initial velocity of the UAV as the positive direction of the vertical axis Y, and the right side of the fuselage as the positive direction of the horizontal axis X, establish a Cartesian coordinate system OXY for the UAV;

[0185] (434) Analyze the reachable space of the trajectory according to different motion stages;

[0186] Based on the coordinate system established in step (433), for the turning stage of the UAV, solve the coordinates of the UAV after the turn according to the direction of the turning angle in step (432); for the acceleration stage of the UAV, it is assumed to accelerate from 10 m / s to 30 m / s at 6 m / s 2 ; for the uniform flight stage of the UAV, fly at a constant speed with the maximum flight speed; solve the position of the UAV at the prediction moment, that is, the final position in the uniform flight stage;

[0187] (435) Discretize the boundary of the reachable space of the trajectory;

[0188] The boundary of the reachable space of the UAV's trajectory at the prediction moment is generated by the set of boundary coordinates of the UAV at that moment under the turning angle in the interval (-360, 360). Discretize the turning angle, sample the turning angle at intervals of 5 degrees, and generate the reachable space of the trajectory with the end point of the previous flight stage as the starting point of the next stage;

[0189] (436) Clear the overlapping area of the reachable space of the trajectory;

[0190] Solve the magnitude of the turning angle of the acceleration trajectory when the first intersection appears in the flight trajectory of the UAV during the turning and acceleration operation time length:

[0191]

[0192] In the formula, s a is the flight distance of the UAV during the acceleration stage, and R min is the minimum turning radius of the UAV; clear the reachable space of the trajectory where overlapping occurs during large-angle turning. The appearance of the overlapping area is due to the absolute value of the turning angle being greater than β0 The set of boundary coordinates X d = {(x c，β , y c，β ), β ∈ (-360, -β 0 ) ∪ (β 0 , 360), the set of boundary coordinates X 0 where the absolute value of the turning angle is less than or equal to β k = {(x c，β , y c，β )}, β ∈ [-β 0 , β 0 , is within the generated closed space. Therefore, it is necessary to remove the set of boundary coordinates when the absolute value of the turning angle is greater than β 0 .

[0193] (437) Regenerate the reachable space of the trajectory for different motion stages;

[0194] Based on the turning angle interval obtained in step (436), sample at intervals of 5 degrees based on this interval, and generate the reachable space of the trajectory covered at different stages when the initial motion prediction duration is 5s and 30s as Figures 3a - 4c . The end point of the previous flight stage is the starting point of the next stage. For a rotary-wing UAV, the moving range in the turning stage is smaller than that in the acceleration and uniform motion stages. As the prediction duration increases, the proportion of the displacement loss of the UAV before reaching the maximum flight speed caused by the turning time compared to the total displacement in the uniform motion stage will be less, that is, the larger the prediction time span, the closer the boundary of the reachable space of the UAV's trajectory is to a circle.

[0195] (5) Index the grid coordinates covered by the trajectory: Approximately represent the reachable space of the trajectory generated in step (4) as the set of grids it covers in the rasterized airspace, and solve the coordinates of the covered grids in the Cartesian coordinate system of the airspace;

[0196] (51) Index the grid coordinates covered by the trajectory in the initial hover state;

[0197] (511) Generate the initial coverage space;

[0198] Generate the initial coverage space according to the circumscribed square of the circle in step (42), and save the coordinate sets of the four vertices of the upper left, lower left, upper right, and lower right of the circumscribed square;

[0199] (512) Expand the initial coverage space;

[0200] Based on the possible incomplete grid coverage at the boundary of the initial coverage space in step (511), round down the grid coordinates of the lower boundaries in the x and y directions, and round up the grid coordinates of the upper boundaries to obtain the expanded coverage space; as Figure 5 shown, the reachable space of the trajectory Sp，t and the initial coverage space S p，t，1 and the minimum integerized grid S p，t，2 are in a relationship of layer-by-layer inclusion;

[0201] (513) Solve the coordinate set of the grid set in the gridified airspace Cartesian coordinate system LL;

[0202] Traverse the grids included in the coverage space in step (512), and retain the grids whose distance from the grid center to the initial position of the UAV is less than or equal to the maximum flight distance, to obtain the coordinate set of the grid set in the gridified airspace Cartesian coordinate system LL;

[0203] (52) Initial motion state trajectory coverage grid coordinate index;

[0204] (521) Construct a discrete boundary coordinate set;

[0205] Based on the motion state of the UAV at the initial moment, obtain the discrete boundary coordinate set of the space reachable by its trajectory at the predicted moment in the UAV Cartesian coordinate system OXY in step (433);

[0206] (522) Coordinate system transformation;

[0207] For each discrete boundary coordinate (x c，β , y c，β ) in the UAV Cartesian coordinate system OXY in step (521), convert it to the coordinate (x L c，β , y L c，β ) in the airspace Cartesian coordinate system LL where the grid is located in step (1):

[0208]

[0209] In the formula, α L，o represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. The initial position coordinates of the UAV are (x o , y o ), to obtain the boundary coordinate set on the airspace Cartesian coordinate system LL;

[0210] (523) Obtain the coverage space;

[0211] Generate the initial coverage space according to step (511), and then expand the initial coverage space according to step (512), that is, obtain the coverage space included in the minimum integerized grid of the space reachable by the trajectory;

[0212] (524) Curve fitting circle equation;

[0213] Based on the trajectory reachable space obtained in step (43), transfer the grid center coordinates (x L ck , y L ck ) to (x ck , y ck ) in the UAV Cartesian coordinate system OXY:

[0214]

[0215] where α L，O represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. The UAV position coordinates at the initial moment are (x o , y o ). Furthermore, calculate the angle β ck between the line connecting the point (x ck , y r，o ) in the UAV Cartesian coordinate system OXY and the initial position O(0, 0) of the UAV and the positive Y-axis:

[0216]

[0217] where (x ck , y ck ) is the coordinate of the center of any grid ck in the UAV Cartesian coordinate system OXY, and β r，o ∈[-180, 180]. Furthermore, construct an arc fitting model fit by the least squares method:

[0218] (x ck，c , y ck，c ) = fit({x c，β , y c，β}), β ∈ [β r，o - ε, β r，o + ε] (9)

[0219] The model input is the set of boundary coordinates X k，r = {(x c，β , y c，β )}, β ∈ [β r，o - ε, β r，o + ε], where ε is used to define the angle range. Consider the points within the small angle range in the set of boundary coordinates X k，r = {(x c，β , y c，β )} as an arc, and perform curve fitting by the least squares method. The output is the center (x ck，c , y ck，c ) of the arc and the radius r ck，c , obtaining the fitting circle equation corresponding to the minimum error;

[0220] (525) Determine the positional relationship between the grid center coordinates and the fitted circle;

[0221] Calculate the Euclidean distance between the grid center coordinates and the center of the fitted circle in step (524). By comparing the magnitude of this Euclidean distance with the radius of the fitted circle, determine whether the grid center coordinates are located inside the fitted circular arc; consider the grids located inside the arc as being within the reachable space of the trajectory.

[0222] (526) Solve for the coordinate set of the grid set in the grid - based airspace Cartesian coordinate system LL;

[0223] Traverse the grids contained in the covered space in step (523), and determine the positional relationship between the grid center coordinates and the fitted circle through step (525) to obtain the coordinate set of the grid set in the grid - based airspace Cartesian coordinate system LL;

[0224] (53) Represent the indexed grid coordinates;

[0225] Represent the coordinate sets obtained by indexing the grid sets in different flight states as the set {ck}, where ck = (i k , j k ) indicates that the grid ck in the grid set is the i - th k grid in the x - direction and the j - th k grid in the y - direction in the airspace Cartesian coordinate system LL.

[0226] (6) Generate the trajectory probability distribution: Solve for the probability that the UAV is located at each grid in the airspace Cartesian coordinate system at the monitoring moment to obtain the corresponding probability distribution of the grid set;

[0227] (61) Solve for the trajectory probability distribution in the initial hovering state;

[0228] (611) Coordinate system translation;

[0229] Consider the subsequent movement of the UAV in the initial hovering state as an undirected Brownian motion, and translate the reachable space of the trajectory from the original (x, y) in the airspace Cartesian coordinate system LL to a new coordinate system (x′, y′) with the position of the UAV at the initial moment as the origin:

[0230]

[0231] where (x 0 , y 0 ) is the position coordinate of the UAV at the initial moment;

[0232] (612) Solve for the trajectory probability distribution;

[0233] In the translated coordinate system after step (611), the predicted UAV position coordinates are obtained through formula (3). The probability density function of k k is obtained. The grid coordinates in the grid set of step (513) are transformed through formula (10). For the grid ck = (i k ), the coordinates of the grid center, lower left, upper left, upper right, and lower right vertices in the translated coordinate system are successively {(x ck ′(t), y ck ′(t)), (x ck，l ′(t), y ck，l ′(t)), (x ck，l ′(t), y ck，u ′(t)), (x ck，r ′(t), y ck，u ′(t)), (x ck，r ′(t), y ck，l ′(t))}. Solve for the position of the UAV at the predicted time, i.e., time t, within the reachable space of the trajectory The probability p 1，ck (t) that it is located within the grid ck:

[0234]

[0235] In the formula, x ck，l ′(t) and x ck，r ′(t) are the abscissas of the upper left and lower left, upper right and lower right vertices of the grid respectively, and y ck，l ′(t) and y ck，u ′(t) are the ordinates of the lower left and lower right, upper left and upper right vertices of the grid respectively. is the probability density function of the coordinates at time t in the initial hovering state ;

[0236] (613) Update the trajectory probability distribution;

[0237] Sum the probability values in each grid obtained in step (612) to get the total probability p 1 (t). p 1 (t) is extremely close to 1 when the grid granularity is fine enough:

[0238]

[0239] In the formula, p 1，ck (t) is the probability that the UAV in the initial hovering state at time t is located within the grid ck. To extend the model to scenarios with a larger grid granularity, p (t) is updated through formula (12): 1，ck (t):

[0240]

[0241] When obtaining the initial hovering state of the UAV, the probability set of the UAV at each grid in the grid set at time t is obtained;

[0242] (62) Solving the probability distribution of the initial motion state trajectory;

[0243] (621) Solving the expectation of the reachable space of the trajectory at the prediction time;

[0244] The subsequent motion of the UAV with the initial motion state is regarded as the combination of undirected Brownian motion and directed Brownian motion. Using the trajectory prediction model based on data migration in step (32), the trajectory prediction value of the UAV at the prediction time is solved, and it is regarded as the expectation E=(μ x (t), μ y (t));

[0245] (622) Translating and rotating the coordinate system;

[0246] Translate the reachable space of the trajectory from the grid-based airspace Cartesian coordinate system LL to the origin with the expected value E of the UAV position at the prediction time obtained in step (621), and rotate it to make the line connecting the initial time coordinates and E the positive X-axis of the new coordinate system EXY. Perform coordinate transformation of the original coordinate system in the new coordinate system:

[0247]

[0248] In the formula, α L，E represents the angle between the airspace Cartesian coordinate system LL and the positive X-axis of the EXY coordinate system;

[0249] (623) Determining the upper and lower limits of the abscissa of the UAV in the reachable space of the trajectory at the prediction time;

[0250] Convert the set of boundary coordinates on the airspace Cartesian coordinate system LL to the coordinate system EXY in step (522), and represent the upper and lower limits of the abscissa of the UAV in the reachable space of the trajectory at the prediction time;

[0251] (624) Determining the upper and lower limits of the ordinate based on the abscissa;

[0252] Calculate the absolute value of the difference between the abscissa in step (623) and the converted set of boundary coordinates, and solve for the index i corresponding to the minimum value in this set x，max and the y corresponding to this index i , if y i ≥0, set the upper limit of the ordinate to y i , execute the loop, set the index i x，max at this position to +∞, and solve for the index i corresponding to the minimum value in this setx，max , until index i x，max corresponding y i < 0, stop the loop, and set the lower limit of the vertical coordinate to y i ; If y i < 0, set the lower limit of the vertical coordinate to y i , execute the loop, set the index i of this set x，max to +∞, and solve for the index i corresponding to the minimum value of this set x，max , until index i x，max corresponding y i ≥0, stop the loop, and set the upper limit of the vertical coordinate to y i ;

[0253] (625) Solve the probability density function at the prediction moment;

[0254] Based on the upper and lower limits of the abscissa determined in step (623) and the upper and lower limits of the ordinate determined in step (624), obtain the probability density function of the UAV position coordinates at the prediction moment when the initial state is in motion through formula (3) within the reachable space of the trajectory;

[0255] (626) Select the integration range of the probability density function;

[0256] Based on the vertex coordinates of the rotated grid in step (622), the binary integration domain will become larger. Calculate the area s of the interior of the true grid integration domain obtained based on the coordinate system rotation angle 1 For example Figure 6a and Figure 6b , the coordinate system rotation angles α L，E are 12° and 50° respectively. When integrating the probability density function established on the EXY coordinate system, selecting the vertex coordinates of the rotated grid ck will cause the binary integration domain to become larger, that is Figure 6a and Figure 6b the area s in 2 , the integration domain of the true grid ck should be Figure 6a and Figure 6b the area s in c，k , establish the relationship between the true grid integration domain s ck and the area s of the binary integration domain of the rotated grid 2 :

[0257]

[0258] (627) Solve the trajectory probability distribution;

[0259] In the translated and rotated coordinate system of step (622), convert the grid coordinates in the grid set of step (526) through formula (14), and the grid ck = (i k , jk ) The coordinates of the grid center, lower left, upper left, upper right, and lower right vertices in the translation coordinate system are successively {(x ck ′(t), y ck ′(t)), (x ck，ll ′(t), y ck，ll ′(t)), (x ck，lu ′(t), y ck，lu ′(t)), (x ck，ru ′(t), y ck，ru ′(t)), (x ck，rl ′(t), y ck，rl ′(t))}. According to the relationships obtained in step (626), establish the double integrals p 1 and p 2 on areas s 2，ck1 (t) and p 2，ck2 (t) respectively:

[0260]

[0261]

[0262] In the formula, (x ck ′(t), y ck ′(t)) is the center coordinate of grid ck, a and b are respectively half of the lengths of two adjacent sides in area s 1 , U ck，x ′(t) and L ck，x ′(t) are respectively the maximum and minimum values of the abscissas of the lower left, upper left, upper right, and lower right vertices of grid ck, and U ck，y ′(t) and L ck，y ′(t)) are respectively the maximum and minimum values of the ordinates of the lower left, upper left, upper right, and lower right vertices of grid ck. is the probability density function of the coordinate at the initial motion state t. According to formula (15), obtain the probability p that the position of the UAV in the reachable space of the trajectory at the initial motion state t is located within grid ck: 2，ck (t):

[0263]

[0264] (628) Solve the probability distribution set;

[0265] Update the probability p 2，ck (t) in step (627) according to the probability update formula (13) in step (613) to obtain the probability set of the UAV located in each grid in the grid set at the initial motion state of the UAV.

[0266] The specific application ways of the present invention are numerous. The above are only the preferred embodiments of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements can still be made, and these improvements should also be regarded as the protection scope of the present invention.

Claims

1. A method for predicting the trajectory distribution of non - cooperative unmanned aerial vehicles (UAVs) based on flight state division, Characterized in that, The steps are as follows: (1) For the selected monitoring airspace range, construct a rasterized airspace; (2) Divide the UAV flight state: According to the flight state of the UAV at the initial moment, divide the flight state of the UAV into a hovering state and a moving state, and set the parameters for predicting the UAV trajectory distribution; (3) Screen the similar trajectory data set and construct a trajectory prediction model based on data migration; (4) Generate the trajectory reachable space: According to the UAV flight state divided in step (2), solve the set of all possible trajectory positions of the UAV at the prediction moment as a region with a closed boundary, and generate the trajectory reachable space of the UAV at the prediction moment; (5) Index the grid coordinates covered by the trajectory: Approximately represent the trajectory reachable space generated in step (4) in the rasterized airspace as the set of grids it covers, and solve the coordinates of the covered grids in the airspace Cartesian coordinate system; (6) Generate the trajectory probability distribution: Solve the probability that the UAV is located in each grid of the airspace Cartesian coordinate system at the monitoring moment to obtain the probability distribution corresponding to the grid set; The specific process of step (6) is as follows: (61) Solve the trajectory probability distribution in the initial hovering state; (611) Translate the coordinate system; (612) Solve the trajectory probability distribution; (613) Update the trajectory probability distribution; Sum the probability values within each grid obtained according to step (612) to obtain the total probability p 1 (t), p 1 (t) is extremely close to 1 when the grid is fine enough: where p 1，ck (t) is the position of the UAV in the reachable space of the trajectory at time t in the initial hovering state the probability of being located within the grid ck. To extend the model to scenarios with a larger grid granularity, p 1，ck (t) is updated by Equation (12): Obtain the probability set that the UAV is located in each grid of the grid set at time t when the UAV is in the initial hovering state; (62) Solve the trajectory probability distribution in the initial moving state; (621) Solve the expectation of the trajectory reachable space at the prediction moment; (622) Translate and rotate the coordinate system; Translate the trajectory reachable space from the rasterized airspace Cartesian coordinate system LL to the origin with the expected value E of the UAV position at the prediction moment obtained in step (621), and rotate it so that the line connecting the coordinates at the initial moment and E is the positive direction of the X - axis of the new coordinate system EXY, and perform the coordinate transformation of the original coordinate system in the new coordinate system: where α L，E represents the angle between the positive direction of the X-axis of the EXY coordinate system and the LL of the spatial Cartesian coordinate system; (623) Determine the upper and lower limits of the abscissa of the UAV in the trajectory reachable space at the prediction moment; (624) Determine the upper and lower limits of the ordinate based on the abscissa; Calculate the absolute value of the difference between the abscissa and the set of transformed boundary coordinates in calculation step (623), and solve for the index i corresponding to the minimum value in this set x，max and the y corresponding to this index i . If y i ≥0, set the upper limit of the ordinate to y i , execute a loop, set the set index i x，max to +∞, solve for the index i corresponding to the minimum value in this set x，max , until the index i x，max corresponding y i <0, stop the loop, and set the lower limit of the ordinate to y i ; if y i <0, set the lower limit of the ordinate to y i , execute a loop, set the set index i x,max to +∞, solve for the index i corresponding to the minimum value in this set x,max , until the index i x,max corresponding y i ≥0, stop the loop, and set the upper limit of the ordinate to y i ; (625) Solve the probability density function at the prediction moment; (626) Select the integration range of the probability density function; Rotating the vertices of the grid according to step (622) will result in an enlarged binary integration domain, and calculate the area s obtained based on the coordinate system rotation angle within the true grid integration domain 1 , and establish the true grid integration domain s ck and the area s of the binary integration domain of the rotated grid 2 The relationship between them is: (627) Solve the trajectory probability distribution; In the translated and rotated coordinate system of step (622), the grid coordinates in the grid set of step (5) are transformed by formula (14). The grid center, lower left, upper left, upper right, and lower right vertices of grid ck = (i k , j k ) in the translated coordinate system are successively {(x ck ′(t), y ck ′(t)), (x ck,ll ′(t), y ck,ll ′(t)), (x ck,lu ′(t), y ck,lu ′(t)), (x ck,ru ′(t), y ck,ru ′(t)), (x ck,rl ′(t), y ck,rl ′(t))}. According to the relationships obtained in step (626), the double integrals p 1 and p 2 on areas s 2,ck1 (t) and s 2,ck2 (t) are respectively established: where, (x ck ′(t), y ck ′(t)) is the center coordinate of grid ck, a and b are respectively half of the two adjacent side lengths of area s 1 , U ck,x′ (t) and L ck,x′ (t) are respectively the maximum and minimum values of the abscissas of the four vertices at the lower left, upper left, upper right, and lower right of grid ck, U ck,y′ (t) and L ck,y′ (t)) are respectively the maximum and minimum values of the ordinates of the four vertices at the lower left, upper left, upper right, and lower right of grid ck, is the probability density function of the coordinate at the initial motion state at time t, and the position of the UAV within the trajectory reachable space at the initial motion state at time t is obtained according to formula (15) The probability p that it is located within grid ck 2,ck (t): (628) Solve the probability distribution set; Update the probability p(t) in step (627) according to the probability update formula (13) in step (613), and obtain the probability set of the UAV located in each grid in the grid set when the UAV is in the initial motion state. 2,ck (t) to obtain the probability set of the UAV located in each grid in the grid set when the UAV is in the initial motion state.

2. The method for predicting the trajectory distribution of non - cooperative UAVs based on flight state division according to claim 1, Characterized in that, In step (1), for the monitoring airspace range, an airspace Cartesian coordinate system LL is established with the lower - left corner as the origin. By setting the number of grids on the x - axis and y - axis, the airspace is divided into a grid set with unique coordinates. The grid coordinates (i, j) indicate that the grid is the i - th grid in the x - axis direction and the j - th grid in the y - axis direction of the airspace. The information stored in each grid includes: the coordinates of the grid center point and the four vertexes in the airspace Cartesian coordinate system LL.

3. The method for predicting the trajectory distribution of non - cooperative UAVs based on flight state division according to claim 2, Characterized in that, The specific process of step (2) is as follows: (21) Divide the UAV flight state; According to the flight state of the UAV at the initial moment, the flight state of the UAV is divided into a hovering state and a motion state; (22) Set the prediction parameters of the UAV trajectory distribution; (221) Set the initial flight speed of the UAV; For the initial speed of the UAV during the trajectory prediction process, set the absolute value of the initial speed in the hovering state to zero, and the absolute value of the initial speed in the motion state to be greater than zero; (222) Set the maximum flight ground speed and the maximum horizontal acceleration of the UAV; According to the factory performance parameters of the drone, it is estimated that the maximum allowable flight speed of the drone under calm wind is less than 30 m / s. The maximum ground speed is set at 30 m / s, and the maximum horizontal acceleration is 6 m / s 2 ; (223) Set the UAV trajectory prediction duration; According to the maximum flight ground speed, maximum horizontal acceleration, maximum static wind speed of the UAV set in step (222) and the pre-time requirement for UAV dangerous behavior recognition, set the prediction duration to be greater than the time required for the UAV to accelerate from the hovering state to the maximum flight ground speed at the maximum horizontal acceleration.

4. According to the non-cooperative UAV trajectory distribution prediction method based on flight state division described in claim 1, it is characterized in that the specific process of the said step (3) is as follows: (31) Screen the similar trajectory dataset; Collect the trajectory segments of non-cooperative UAVs. According to the existing cooperative UAV trajectory dataset, express the potential characteristics of the UAV trajectory through the flight speed variance and the cumulative heading change amount, and further describe the stability of the flight state of each trajectory segment. Based on the similarity of the flight speed variance and cumulative heading change amount values, measure the distance between trajectories, and screen out the trajectory dataset similar to the non-cooperative UAV trajectory; (32) Construct a trajectory prediction model based on data migration; Use the similar trajectory dataset screened in step (31) as the training sample of the D-GRU trajectory prediction model, and regard the trained D-GRU trajectory prediction model as the trajectory prediction model based on data migration.

5. According to the non-cooperative UAV trajectory distribution prediction method based on flight state division described in claim 4, it is characterized in that the specific process of the said step (4) is as follows: (41) Space-time constrained Brownian bridge; Regarding the motion characteristics of the UAV in different flight states divided in step (21), when representing the position of the UAV at time t as (x(t), y(t)) according to the truncated Brownian bridge, x(t) follows a truncated probability density function: At time t, the UAV reaches the coordinate within the trajectory reachable space For any y(t) follows a truncated probability density function: In addition, the probability that the drone reaches the coordinate within the reachable space of the trajectory is as follows: Wherein, (x(t), y(t)) is the position of the UAV at time t. For the normal distribution X ~ N(μ, σ 2 ), is the probability density function, Φ(t) is the distribution function, U x (t), U y (t) are the upper boundaries in the x and y directions, L x (t), L y (t) are the lower boundaries in the x and y directions; (42) Trajectory reachable space of the UAV in the initial hovering state; For the UAV in the hovering state at the initial moment, the trajectory reachable space at the prediction moment is a circle, with the center at the position at the initial moment and the radius being the maximum flight distance to the prediction moment. The maximum flight distance includes two parts: the accelerated flight distance and the distance of flying at the maximum flight speed at a constant speed; (43) Trajectory reachable space of the UAV in the initial motion state; (431) Determine the minimum turning radius; Obtain its minimum turning radius according to the maximum overload coefficient allowed in the design of the UAV: where R min is the minimum turning radius of the UAV, g represents the acceleration due to gravity, |v 0 | represents the magnitude of the flight speed of the UAV at the start of turning, n max is the maximum overload coefficient. It is assumed that the UAV with an initial velocity starts turning at the initial moment, and during the turning, the initial velocity remains constant in magnitude and only changes in direction; (432) Set the direction of the turning angle; Define the difference between the UAV speed directions before and after turning as the turning angle. When turning right, the turning angle is positive; when turning left, the turning angle is negative; when keeping the initial speed direction unchanged, the turning angle is 0; (433) Establish the UAV Cartesian coordinate system; Taking the initial position of the UAV as the coordinate origin, the direction of the initial velocity of the UAV as the positive direction of the vertical axis Y, and the right side of the fuselage as the positive direction of the horizontal axis X, a Cartesian coordinate system OXY for the UAV is established; (434) Analyze the reachable space of the trajectory according to different motion stages; Based on the coordinate system established in step (433), for the turning stage of the UAV, solve the coordinates of the UAV after the turn according to the direction of the turning angle in step (432); for the acceleration stage of the UAV, it is set to accelerate from the initial velocity to the maximum flight speed with the maximum horizontal acceleration; for the uniform flight stage of the UAV, it flies at a constant speed with the maximum flight speed; solve the position of the UAV at the prediction moment, that is, the last position in the uniform flight stage; (435) Discretize the boundary of the reachable space of the trajectory; The boundary of the reachable space of the UAV's trajectory at the prediction moment is generated by the set of boundary coordinates of the UAV at this moment when the turning angle is in the interval (-360, 360). Discretize the turning angle, sample the turning angle at equal interval degrees, and generate the reachable space of the trajectory with the end point of the previous flight stage as the starting point of the next stage; (436) Clear the overlapping area of the reachable space of the trajectory; Solve the magnitude of the turning angle of the acceleration trajectory when the first intersection appears in the flight trajectory of the UAV during the time length of turning and accelerating operation: where s a is the flight distance during the acceleration phase of the UAV, and R min is the minimum turning radius of the UAV; clear the accessible space of the trajectory that overlaps during turning at too large an angle; (437) Regenerate the reachable space of the trajectory in different motion stages; Sample the turning angle interval obtained in step (436) at equal interval degrees to generate the reachable space of the trajectory covered by the initial motion in different stages.

6. The non-cooperative UAV trajectory distribution prediction method based on flight state division according to claim 5, characterized in that the specific process of step (5) is as follows: (51) Index the grid coordinates covered by the trajectory in the initial hover state; (511) Generate the initial coverage space; Generate the initial coverage space according to the circumscribed square of the circle in step (42), and save the coordinate sets of the upper left, lower left, upper right, and lower right vertices of the circumscribed square; (512) Expand the initial coverage space; Based on the possible incomplete grid coverage at the boundary of the initial coverage space in step (511), round down the grid coordinates of the lower boundary in the x and y directions and round up the grid coordinates of the upper boundary to obtain the expanded coverage space; (513) Solve the coordinate set of the grid set in the grid-based airspace Cartesian coordinate system LL; Traverse the grids included in the coverage space in step (512), and retain the grids whose distance from the grid center to the initial position of the UAV is less than or equal to the maximum flight distance to obtain the coordinate set of the grid set in the grid-based airspace Cartesian coordinate system LL; (52) Index the grid coordinates covered by the trajectory in the initial motion state; (521) Construct a discrete boundary coordinate set; Based on the motion state of the UAV at the initial moment, obtain the discrete boundary coordinate set of the reachable space of its trajectory at the prediction moment in the UAV Cartesian coordinate system OXY in step (433); (522) Coordinate system transformation; For each discrete boundary coordinate (x c，β , y c，β ) in the UAV Cartesian coordinate system OXY of step (521), convert it to the coordinate (x L c，β , y L c，β ) in the airspace Cartesian coordinate system LL where the grid is located in step (1): where α L，o represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. At the initial moment, the position coordinates of the UAV are (x o , y o ), and a set of boundary coordinates on the airspace Cartesian coordinate system LL is obtained; (523) Obtain the coverage space; Generate the initial coverage space according to step (511), and then expand the initial coverage space according to step (512), that is, obtain the coverage space included in the minimum integerized grid of the trajectory reachable space; (524) Curve fitting circle equation; Based on the reachable space of the trajectory obtained in step (43), transfer the grid center coordinates (x L ck , y L ck ) to (x ck , y ck ) in the UAV Cartesian coordinate system OXY: Where α L,O represents the angle of counterclockwise rotation from the airspace Cartesian coordinate system LL to the UAV Cartesian coordinate system OXY. At the initial moment, the UAV position coordinates are (x o , y o ). Then, calculate the angle β ck , y ck ) in the UAV Cartesian coordinate system OXY between the line connecting the point (x r,o : Where, (x ck , y ck ) is the coordinate of the center of any grid ck in the Cartesian coordinate system OXY of the UAV, β r,o ∈ [-180, 180], and then a circular arc fitting model fit is constructed by the least squares method: (x ck,c ,y ck,c ) = fit({x c,β ,y c,β}), β ∈ [β r,o - ε, β r,o + ε] (9) The model input is the set of boundary coordinates X at the prediction moment of the UAV k,r ={(x c,β ,y c,β ),}, β ∈ [β r,o -ε, β r,o +ε], ε is used to define the angular range. The points within the small angular range in the set of boundary coordinates X k,r ={(x c,β ,y c,β )} are regarded as an arc. Curve fitting is performed by the least squares method, and the output is the center (x ck,c , y ck,c ) of the arc and the radius r ck,c , obtaining the fitting circle equation corresponding to the minimum error; (525) Determine the positional relationship between the grid center coordinates and the fitted circle; Calculate the Euclidean distance between the grid center coordinates and the center of the fitted circle in step (524). By comparing the size of this Euclidean distance with the radius of the fitted circle, determine whether the grid center coordinates are located inside the fitted circular arc; regard the grids located inside the circular arc as being inside the trajectory reachable space; (526) Solve the coordinate set of the grid set in the grid-based airspace Cartesian coordinate system LL; Traverse the grids contained in the coverage space in step (523), and determine the positional relationship between the grid center coordinates and the fitted circle through step (525) to obtain the coordinate set of the grid set in the grid-based airspace Cartesian coordinate system LL; (53) Represent the grid coordinates of the index; The coordinate sets obtained by indexing the grid sets in different flight states are represented as the set {ck}, where ck = (i k , j k ) indicates that the grid ck in the grid set is the i-th k grid in the x direction and the j-th k grid in the y direction in the airspace Cartesian coordinate system LL.

7. The non-cooperative UAV trajectory distribution prediction method based on flight state division according to claim 6, characterized in that, the specific process of step (6) is as follows: (61) Solve the trajectory probability distribution in the initial hovering state; (611) Coordinate system translation; The subsequent movement of the UAV in the initial hovering state is regarded as an undirected Brownian motion, and the trajectory reachable space is translated from (x, y) in the original airspace Cartesian coordinate system LL to a new coordinate system (x′, y′) with the UAV position at the initial moment as the origin: where $(x 0 , y 0 )$ is the position coordinate of the UAV at the initial moment; (612) Solve the trajectory probability distribution; In the translated coordinate system of step (611), the position coordinates of the UAV at the prediction moment are obtained through formula (3). The probability density function of, the grid coordinates in the grid set of step (513) are transformed through formula (10), and the grid center of the grid ck = (i k , j k ) and the coordinates of the lower left, upper left, upper right, and lower right vertices in the translated coordinate system are successively {(x ck ′(t), y ck ′(t)), (x ck,l ′(t), y ck,l ′(t)), (x ck,l ′(t), y ck,u ′(t)), (x ck,r ′(t), y ck,u ′(t)), (x ck,r ′(t), y ck,l ′(t))}, and solve for the position of the UAV in the reachable space of the trajectory at the prediction moment, i.e., at moment t The probability p 1，ck (t) that it is located within the grid ck: where x ck,l ′(t) and x ck,r ′(t) are the abscissas of the upper-left and lower-left, upper-right and lower-right vertices of the grid respectively, and y ck,l ′(t) and y ck,u ′(t) are the ordinates of the lower-left and lower-right, upper-left and upper-right vertices of the grid respectively, is the probability density function of the coordinates at time t in the initial hovering state ; (613) Update the trajectory probability distribution; (62) Solve the trajectory probability distribution in the initial motion state; (621) Solve the expectation of the trajectory reachable space at the prediction moment; The subsequent motion of the UAV in the initial motion state is regarded as a combination of undirected Brownian motion and directed Brownian motion. The trajectory prediction value of the UAV at the prediction time is solved by using the trajectory prediction model based on data migration in step (32), and it is regarded as the expectation E=(μ x (t), μ y (t)) of the trajectory reachable space at the prediction time; (622) Translate and rotate the coordinate system; (623) Determine the upper and lower limits of the abscissa of the UAV in the trajectory reachable space at the prediction moment; Convert the boundary coordinate set on the airspace Cartesian coordinate system LL to the coordinate system EXY in step (522) to represent the upper and lower limits of the abscissa of the UAV in the trajectory reachable space at the prediction moment; (624) Determine the upper and lower limits of the ordinate based on the abscissa; (625) Solve the probability density function at the prediction moment; According to the upper and lower limits of the abscissa determined in step (623) and the upper and lower limits of the ordinate determined in step (624), the predicted position coordinates of the drone at the prediction moment in the initial moving state are obtained by formula (3). The probability density function within the trajectory reachable space; (626) Select the integration range of the probability density function; (627) Solve the trajectory probability distribution; (628) Solve the probability distribution set.

Images