A method for conflict risk assessment of non - cooperative unmanned aerial vehicles invading cooperative unmanned aerial vehicles

Abstract

The present invention discloses a method for conflict risk assessment of non - cooperative drones invading cooperative drones, including: collecting real - trajectory data of various types of drones to generate trajectory time - series slices; identifying conflicts according to the predicted information of non - cooperative drones and cooperative drones at a certain future moment and the set protected area range; performing simulation on the conflict - existing point - pair combinations to generate Monte Carlo simulation samples; evaluating the possibility and severity of the occurrence of dangerous behaviors of non - cooperative drones invading cooperative drones, determining risk assessment indicators and rating them; and conducting conflict risk assessment and level calculation for the scenario of non - cooperative drones invading cooperative drones. The method of the present invention can conduct risk assessment on the conflict scenario between non - cooperative drones and a single drone and quantify the conflict risk level, thereby providing more accurate alarm information.

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 evaluating the conflict risk of non-cooperative UAVs invading cooperative UAVs. Background Art

[0002] In actual UAV control tasks, simply achieving real-time UAV trajectory prediction is difficult to directly assist in the selection of conflict resolution strategies. Therefore, on the basis of UAV trajectory prediction, identifying conflicts between UAVs and evaluating the conflict risk has become an important research direction in UAV safety management.

[0003] Regarding the risk level classification of light and smaller UAVs in the low-altitude airspace, there is currently no unified and standardized regulation and development standard at home and abroad, and most of the UAVs in the airspace are recreational UAVs flying freely or autonomous operation UAVs; relatively speaking, they may have characteristics such as non-cooperativeness and dynamic uncertainty, which makes it difficult to detect the conflict risk of low-altitude UAVs, and at the same time affects the reliability and efficiency of warning of dangerous behaviors. Therefore, for the scenario of conflicts between non-cooperative UAVs and cooperative UAVs, based on the motion characteristics of UAVs and the protected area range set by cooperative UAVs, identify conflicts, evaluate and quantify the conflict risk level, achieve a reasonable classification of different conflict risk levels, provide alarm information in a timely manner, and avoid the occurrence of conflict behaviors. 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 evaluating the conflict risk of non-cooperative UAVs invading cooperative UAVs to solve the problems in the existing technologies that it is difficult to evaluate the conflict risk between non-cooperative UAVs and cooperative UAVs and the risk level is not clear; the method of the present invention can evaluate the risk of conflicts between non-cooperative UAVs and a single UAV and quantify the conflict risk level, so as to provide more accurate alarm information.

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

[0006] A method for evaluating the conflict risk of non-cooperative UAVs invading cooperative UAVs according to the present invention comprises the following steps:

[0007] 1) Collect real trajectory data of various types of UAVs, and perform trajectory data processing and segmentation to generate trajectory time series slices;

[0008] 2) Identify conflicts according to the prediction information of non-cooperative UAVs and cooperative UAVs at a certain future moment and the set protected area range;

[0009] 3) Based on the conflict identification result in step 2), perform simulation on the conflicting point pair combinations to generate Monte Carlo simulation samples;

[0010] 4) Based on the simulation samples generated in step 3), evaluate the possibility and severity of the occurrence of dangerous behaviors of non-cooperative drones invading cooperative drones, determine risk assessment indicators and conduct ratings;

[0011] 5) Based on the risk assessment indicators determined in step 4), conduct conflict risk assessment and level calculation for the scenario of non-cooperative drones invading cooperative drones.

[0012] Furthermore, the specific process of step 1) is as follows:

[0013] Collect real trajectory data of various types of drones. Each trajectory point in the trajectory data includes timestamp, identity ID, latitude, longitude, flight relative height, speed, and attitude angle information. Convert the longitude, latitude, and height in the WGS-84 coordinate system into (x, y, z) in the Cartesian coordinate system. Perform equidistant sampling on the real trajectory data of drones, use the sliding time window method to slice and divide the trajectory sequence data of drones, and perform normalization processing to generate trajectory time series slices of various types of drones.

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

[0015] 21) Judge the relative position relationship between the predicted position of a non-cooperative drone at a future moment and the conflict protection area of a cooperative drone;

[0016] 211) Define the conflict protection area and collision protection area of the cooperative drone;

[0017] Taking the predicted position of the cooperative drone at a future moment as the midpoint of the horizontal and vertical directions of the protection area, considering the two processes of conflict and collision, set the shapes of the conflict protection area and the collision protection area as cylinders, and the conflict protection area is formed by expanding the collision protection area outward; the size calculation formula of the collision protection area is as follows:

[0018]

[0019]

[0020] In the formula, h is the height of the cylinder of the collision protection area, r is the radius of the cylinder of the collision protection area, h u is the height of the external dimension of the drone, r u is half of the maximum value of the horizontal dimension of the drone, G y and G h are the vertical and horizontal hovering precisions of the drone in the GPS mode respectively; the cylinder height D ver and radius D horEnlarge it by a certain proportion based on the height and radius of the cylinder in the collision protection area;

[0021] 212) Determine the conflict interval of the UAV;

[0022] Determine whether the predicted three-dimensional position coordinates of UAV B as a particle at the t-th moment are both within the collision protection area of cooperative UAV A in the horizontal and vertical directions. If so, the following formulas (3) and (4) hold simultaneously, and a conflict occurs between UAV A and UAV B. The specific judgment formulas are as follows:

[0023]

[0024] |P(z) i-coo -P(z) i-non |<D ver (4)

[0025] In the formula, P(x) i-coo 、P(y) i-coo 、P(z) i-coo are the x coordinate in the horizontal direction, the y coordinate in the horizontal direction, and the z coordinate in the vertical direction of the cooperative UAV respectively, and P(x) i-non 、P(y) i-non 、P(z) i-non are the x coordinate in the horizontal direction, the y coordinate in the horizontal direction, and the z coordinate in the vertical direction of the non-cooperative UAV respectively;

[0026] If both formulas (3) and (4) do not hold, then there is no conflict between UAV A and UAV B; if only one of formulas (3) and (4) holds, then go to step 22);

[0027] 22) Based on the judgment result in step 21), if formulas (3) and (4) are not satisfied simultaneously, conflict identification is further carried out in combination with the predicted flight intention of the non-cooperative UAV at future moments; among them, the flight intention within the time period from a certain future moment to the next moment is characterized by three aspects, including horizontal flight behavior, vertical flight behavior, and acceleration situation; the horizontal flight behavior is specifically the course angle interval between two adjacent moments, which is divided into four intervals at intervals of 90°; the vertical flight behavior is specifically the pitch angle interval between two adjacent moments, including three intervals of descending, climbing, and level flight; the acceleration situation consists of three situations of acceleration, uniform speed, and deceleration between two adjacent moments; if the predicted value of the three-dimensional position coordinates of the non-cooperative UAV B at the t-th moment is only within the conflict protection area of the cooperative UAV A in the horizontal or vertical direction, and the distance between the two in the other direction is less than the flight distance traveled at the maximum speed under a single step of the UAV, then it is judged whether the probability that the flight behavior interval of the non-cooperative UAV B in the corresponding direction is in the same interval as the cooperative UAV A exceeds a predetermined threshold. If it exceeds, it means that there is a conflict risk between the predicted trajectories of the non-cooperative UAV B and the cooperative UAV A at the t-th moment; if it does not exceed, there is no conflict between the cooperative UAV A and the non-cooperative UAV B.

[0028] 23) Based on the results in step 21) and step 22), obtain the prediction label of the conflict behavior between the non-cooperative UAV and the cooperative UAV at the t-th moment; meanwhile, when the relative distances in the horizontal and vertical directions between the true position of the non-cooperative UAV at the t-th moment and the true trajectory coordinates of the cooperative UAV are less than or equal to the set conflict protection area distance, the true label of the conflict event is yes, otherwise it is no.

[0029] 24) Conduct an effectiveness evaluation of the conflict identification result.

[0030] 241) Calculate the conflict identification accuracy.

[0031] The conflict identification accuracy is the proportion Acc of the number of correctly predicted samples to the total number of samples. The specific calculation formula is as follows:

[0032]

[0033] In the formula, y l , y' l are the true label and the predicted label of sample l respectively, and n is the number of predicted samples.

[0034] 242) Calculate the precision and recall of the conflict identification.

[0035] In the conflict identification of the non-cooperative UAV intrusion scenario, the specific calculation formulas for precision and recall are as follows:

[0036]

[0037]

[0038] Wherein, R is the recognition recall rate, and P is the recognition precision rate;

[0039] 243) Calculation of the comprehensive evaluation index for conflict recognition;

[0040] The emphasis degree of the recognition task on the precision rate P and the recall rate R is characterized by the balance factor α. When α > 1, it emphasizes the recall rate R; when α < 1, it emphasizes the precision rate P; P and R are synthesized into an evaluation index F α , and the specific calculation formula is as follows:

[0041]

[0042] Set the threshold value of the comprehensive evaluation index as F α * , when F α > F α * , enter step 3); otherwise, return to step 1).

[0043] Furthermore, the specific process of the said step 3) is as follows:

[0044] 31) Establish the kinematic equation of the unmanned aerial vehicle;

[0045] Establish the kinematic equation of the unmanned aerial vehicle to ensure the physical feasibility and authenticity of the randomly generated motion trajectories within each adjacent moment. The specific formula is as follows:

[0046]

[0047] Wherein, a is the acceleration vector, V is the velocity vector, T is the thrust vector, K is the empirical drag coefficient, m is the mass of the unmanned aerial vehicle, and g is the acceleration; among them, K is related to the change of the attitude of the unmanned aerial vehicle. At different inclination angles , the relationship between K and the axial drag coefficient K a , and the radial drag coefficient K s is:

[0048]

[0049] Set the drag term of the unmanned aerial vehicle to be axisymmetric, and solve it by combining the relevant performance parameters of the unmanned aerial vehicle in the normal gear and the sports gear to determine the initial thrust and drag. The specific formula is as follows:

[0050]

[0051] Wherein, and are the maximum tilt angles of the drone in the motion mode and the normal mode respectively, and V1 and V2 are the maximum horizontal speeds of the drone in the motion mode and the normal mode respectively;

[0052] 32) The thrust vector and the attitude inclination angle required for generating a specified velocity vector;

[0053] It is set that at the initial stage of each time step, when the attitude and speed change, the thrust vector of the drone is in the same direction as the attitude vector, and the direction of the drag vector is opposite to the direction of the velocity vector, and the unit attitude vector and the unit velocity vector The specific calculation formula is as follows:

[0054]

[0055]

[0056] In the formula, ψ and θ are the heading angle and the pitch angle respectively, and are the rotation angles around the x-axis and the y-axis (the traveling direction of the drone) respectively;

[0057] The randomly generated speed and attitude of the drone meet the physical feasibility requirements of the motion trajectory, and the specific formula is as follows:

[0058]

[0059] In the formula, D is the drag vector, based on the constraint condition In the case of the acceleration vector a = 0, combining Equation (9) and Equation (14), by iteratively changing and until the thrust vector T is perpendicular to the unit attitude vector the thrust vector and the attitude inclination angle required for generating a given velocity vector are generated, so that the randomly generated speed and attitude meet the physical feasibility requirements of the motion trajectory;

[0060] 33) Generate Monte Carlo simulation samples;

[0061] For each Monte Carlo simulation sample, the flight trajectory randomly generates speed, heading angle and elevation angle at each time step. The initial generation of speed, heading angle and elevation angle follows a normal distribution, and Monte Carlo simulation samples of the drone at different moments are generated respectively. The specific generation method is as follows:

[0062]

[0063] In the formula, |V| represents the magnitude of the speed; |V target | represents the magnitude of the average speed; |V max| represents the magnitude of the maximum speed; ψ and θ are the heading angle and pitch angle respectively; V target , ψ target , θ target are the mean values of speed, heading angle, and pitch angle respectively, and σ V , σ ψ , σ θ are the variances of speed, heading angle, and pitch angle respectively.

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

[0065] 41) Evaluate the likelihood of a non - cooperative UAV invading a cooperative UAV's dangerous behavior;

[0066] 411) Calculate the predicted invasion probability of the UAV;

[0067] The predicted invasion probability represents the probability of a non - cooperative UAV invading the conflict protection area and collision protection area of a cooperative UAV at the prediction moment, and is obtained by calculating the percentage of Monte Carlo simulation samples located within the conflict protection area or collision protection area of the cooperative UAV. The specific calculation formula is as follows:

[0068]

[0069] In the formula, μ1 is the predicted invasion probability, f is the number of Monte Carlo simulation samples, s k ∈I indicates that the position of the k - th Monte Carlo simulation sample generated is within the protection area I of the cooperative UAV c , where c = 1 is the conflict protection area and c = 2 is the collision protection area;

[0070] 412) Calculate the predicted conflict time of the UAV;

[0071] When setting the prediction step size as i, the predicted position of the cooperative UAV is (x i , y i , z i ), and the distribution position of the Monte Carlo simulation samples of the non - cooperative UAV is where, j is the number of generated simulation samples; by calculating the distances from the cooperative UAV to each sample point of the non - cooperative UAV, the sample distance set is obtained, where, is the distance between the predicted point of the cooperative UAV and the p - th simulation sample of the non - cooperative UAV at the prediction step size i. The specific calculation formula is as follows:

[0072]

[0073] The minimum distance between the cooperative UAV and the non - cooperative UAV at the prediction step size i The specific calculation formula is as follows:

[0074]

[0075] It is assumed that the cooperative UAV is relatively stationary and does not perform any maneuvers. The calculation of the predicted conflict time is carried out based on the worst-intention assumption that the non-cooperative UAV heads straight towards the cooperative UAV at its maximum flight speed. The specific calculation formula is as follows:

[0076]

[0077] In the formula, μ2 is the predicted conflict time, is the minimum distance between the cooperative UAV and the non-cooperative UAV at the prediction step i, is the maximum flight speed of the non-cooperative UAV;

[0078] 413) Calculate the predicted conflict intention of the UAV;

[0079] Set the horizontal intention interval of the non-cooperative UAV in the time period from prediction step i to i + 1 as I h , the vertical intention interval I v , and the speed change situation I a When accelerating, I a = 1, when moving at a constant speed, I a = 0, when decelerating, I a = -1; the predicted position of the cooperative UAV at step i is (x i , y i , z i ); judge the conflict intention level by combining the relative position of the intention interval and the speed change situation;

[0080] When (x i , y i ) ∈ I h , z i ∈ I v and I a = 1, the non-cooperative UAV will enter the position interval of the cooperative UAV in the horizontal and vertical directions at the next moment, and when accelerating, it means that the non-cooperative UAV has an intrusion intention and the risk level of an accident is the highest; when and I a = -1, the position intervals of the non-cooperative UAV in the horizontal and vertical directions are both separated from the cooperative UAV at the next moment, and when decelerating, it means that the non-cooperative UAV maintains a safe interval between UAVs and the risk level of an accident is the lowest;

[0081] 414) Calculate the UAV distance avoidance index;

[0082] The position and movement direction of the conflicting entities are used to calculate the drone distance avoidance index. The specific calculation formula is as follows:

[0083]

[0084] Where μ4 is the distance avoidance index, X is the distance between the cooperative UAV and the non-cooperative UAV at the prediction moment, and D mar The minimum distance required for the cooperative UAV to safely perform an evasive maneuver, indicating the minimum avoidance range;

[0085] 415) Calculate the UAV direction avoidance index;

[0086] The geometric definition of the speed barrier method is used to determine whether the cooperative UAV and the non-cooperative UAV are on a potential collision route. A circle is drawn with the position of the non-cooperative UAV at the prediction moment as the center and the horizontal safety interval of the conflict protection zone as the radius. The two tangent segments of the circle are drawn with the position of the cooperative UAV at the prediction moment, and the cone area formed is the collision cone. If the speed V of the cooperative UAV is h The speed V of the non-cooperative drone r The vector sum V s If it is within the collision cone, there is a possibility of conflict. The azimuth avoidance index μ5 is calculated as follows:

[0087]

[0088] Where ω is half of the detection angle of the cooperative UAV, V s represents the velocity vector sum of the cooperative UAV and the non-cooperative UAV, α1 and α2 are V s The angle formed by the two boundaries of the host detection angle;

[0089] 416) Determine the probability of non-cooperative UAV invading the conflict area, the probability of non-cooperative UAV invading the collision area, the expected conflict time between non-cooperative UAV and cooperative UAV, the expected conflict intention of non-cooperative UAV, the distance avoidance index of cooperative UAV, the horizontal direction avoidance index of cooperative UAV, and the vertical direction avoidance index of cooperative UAV as evaluation indicators of the possibility of dangerous behavior, and divide them into five levels from 1 to 5 respectively;

[0090] 42) Assess the severity of dangerous behaviors of non-cooperative drones invading cooperative drones;

[0091] 421) Get the drone crash area;

[0092] The drone crash area is divided into vertical collision area and horizontal collision area. The specific calculation formula of the vertical collision area is as follows:

[0093] Re v =π(r p+r uav ) 2 (22)

[0094] In the formula, Re v is the vertical collision area, r p is the average radius of the human body, and r uav is the maximum size radius of the UAV;

[0095] The calculation formula for the horizontal collision area is as follows:

[0096] Re c = 2(r p + r uav )·d + π(r p + r uav ) 2 (23)

[0097] In the formula, Re c is the horizontal collision area, d represents the horizontal distance traveled by the falling UAV when it descends to the height of the human body, and d = H p / tanγ; H p is the average height of the human body, and γ is the angle formed by the velocity vector and the horizontal plane or the human body surface impact;

[0098] 422) Estimate the ground collision kinetic energy of the UAV;

[0099] Calculate the UAV collision point velocity V imp , obtained through the maximum flight speed V x and the free fall speed V y relative to the maximum flight height, as follows:

[0100]

[0101] Through V x and V y , solve the UAV collision angle γ, as follows:

[0102]

[0103] According to the velocity magnitude V imp of the UAV at the impact moment and the maximum takeoff mass M, the kinetic energy E c of the crashed UAV at the impact point can be solved, as follows:

[0104]

[0105] 423) Calculate the probability of death in a UAV ground collision;

[0106] Consider the probability of death P f under the shielding factor. The specific calculation formula is as follows:

[0107]

[0108] In the formula, P s represents the shielding factor, and E h represents the impact energy corresponding to a 50% probability of death when the shielding factor is 6, and E d represents the lowest impact energy threshold that causes death, and E c is the kinetic energy of the crashed drone at the impact point, and μ represents a correction factor used to improve the estimation of low kinetic energy values less than or near the E d threshold, as follows:

[0109]

[0110] 424) Calculate the number of ground casualties per flight hour;

[0111] Based on formulas (22) - (28), the specific calculation formula for the number of ground casualties per flight hour of the drone is as follows:

[0112] N p = P b ·Re c ·D p ·P f (29)

[0113] In the formula, N P is the number of ground casualties per flight hour, P b is the failure probability of the drone. Since the situation considered is that conflicts or collision events in the intrusion scenario may cause accidents and failures of the drone, therefore, P b is set to μ1, Re c is the horizontal collision area, D p is the population density of the drone operation area, and P f is the probability of death;

[0114] 425) Determine the probability of death and the number of ground casualties per flight hour considering the shielding factor as evaluation indicators for the severity of dangerous behavior occurrences, and divide the probability of death into five levels from 1 to 5.

[0115] Furthermore, the drone distance avoidance index refers to the ability of a cooperative drone to avoid potential conflicts or collision behaviors by performing maneuvers at the relative distance level between the cooperative drone and the non - cooperative drone in a potential conflict scenario.

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

[0117] 51) Set the performance parameters and environmental parameters of the unmanned aerial vehicle (UAV), including the maximum takeoff mass of the UAV, maximum operating speed, maximum size radius, detection distance, detection angle, population density in the commercial area, shielding coefficient in the commercial area, average height of standing human body, average radius of human body, and gravitational acceleration;

[0118] 52) Calculate the risk and determine the warning level for the conflict and collision prediction simulation example cases at prediction steps of 5s, 10s, 15s, and 20s respectively, including the probability of invading the conflict and collision protection area, predicted conflict time, predicted conflict intention, distance avoidance index, azimuth avoidance index in the horizontal and vertical directions, probability of death, and number of ground casualties per flight hour.

[0119] Advantages of the present invention:

[0120] Based on the scenario of non - cooperative UAVs invading cooperative UAVs, the present invention combines the predicted points and kinematic modeling method to generate Monte Carlo simulation trajectory samples, and determines the risk assessment indicators from two angles of the possibility and severity of UAV conflicts, so as to calculate the conflict risk assessment and level estimation results of the conflict scenario. The method of the present invention provides a reference basis for the quantification of the UAV operation alarm level. Description of the Drawings

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

[0122] Figure 2 It is a schematic diagram of the risk level and index parameter values under different prediction steps of the embodiment of the present invention. Detailed Embodiment

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

[0124] Refer to Figure 1 As shown, a method for conflict risk assessment of non - cooperative UAVs invading cooperative UAVs of the present invention is as follows:

[0125] 1) Collect real - trajectory data of various types of UAVs, and perform trajectory data processing and segmentation to generate trajectory time - series slices;

[0126] Collect real - trajectory data of multiple types of drones. Each trajectory point in the trajectory data includes timestamp, identity ID, latitude, longitude, relative flight altitude, speed, and attitude angle information. Convert the longitude, latitude, and altitude in the WGS - 84 coordinate system into (x, y, z) in the Cartesian coordinate system. Perform equally - spaced sampling on the real - trajectory data of drones, use the sliding time - window method to slice and divide the trajectory sequence data of drones, and perform normalization processing to generate trajectory time - series slices of various types of drones.

[0127] In the example, equally - spaced sampling is performed on the drone trajectory data at an interval of 1 s; for the divided trajectory slices, set the first trajectory point as the coordinate origin (0, 0, 0) of the slice, and convert the remaining trajectory points into relative coordinates relative to the first point. Use the displacement sequence relative to the initial position of the slice as the input sample distance feature.

[0128] 2) Identify conflicts based on the predicted information of non - cooperative drones and cooperative drones at a future moment, and the set protected area range; specifically as follows:

[0129] 21) Judge the relative position relationship between the predicted position of the non - cooperative drone and the conflict protected area of the cooperative drone at a future moment.

[0130] 211) Define the conflict protected area and the collision protected area of the cooperative drone.

[0131] Take the predicted position of the cooperative drone at a future moment as the mid - point in the horizontal and vertical directions of the protected area. Considering the two processes of conflict and collision, set the shapes of the conflict protected area and the collision protected area as cylinders. The conflict protected area is formed by expanding the collision protected area outward; the size calculation formula of the collision protected area is as follows:

[0132]

[0133]

[0134] In the formula, h is the height of the cylinder of the collision protected area, r is the radius of the cylinder of the collision protected area, h u is the height of the outer dimension of the drone, r u is half of the maximum value of the horizontal dimension of the drone, G y and G h are the vertical and horizontal hovering precisions of the drone in the GPS mode respectively; the height D ver and radius D hor of the cylinder of the conflict protected area are expanded by a certain proportion based on the height and radius of the cylinder of the collision protected area.

[0135] In the example, set the height D ver and radius Dhor Expand by a factor of 8 based on the height and radius of the cylinder in the collision protection area.

[0136] 212) Determine the conflict interval of the UAV.

[0137] Judge whether the predicted values of the three-dimensional position coordinates of UAV B as a particle at the t-th moment are both within the collision protection area of cooperative UAV A in the horizontal and vertical directions. If so, the following formulas (3) and (4) hold simultaneously, and a conflict occurs between UAV A and UAV B. The specific judgment formulas are as follows:

[0138]

[0139] |P(z) i-coo -P(z) i-non |<D ver (4)

[0140] In the formula, P(x) i-coo 、P(y) i-coo 、P(z) i-coo are the x coordinate in the horizontal direction, the y coordinate in the horizontal direction, and the z coordinate in the vertical direction of the cooperative UAV respectively, and P(x) i-non 、P(y) i-non 、P(z) i-non are the x coordinate in the horizontal direction, the y coordinate in the horizontal direction, and the z coordinate in the vertical direction of the non-cooperative UAV respectively;

[0141] If both formulas (3) and (4) do not hold, then there is no conflict between UAV A and UAV B; if only one of formulas (3) and (4) holds, then go to step 22);

[0142] 22) Based on the judgment result of step 21), if formulas (3) and (4) are not satisfied simultaneously, conflict identification is further carried out in combination with the predicted flight intention of the non-cooperative UAV at the future moment; among them, the flight intention within the time period from a certain future moment to the next moment is characterized by three aspects, including horizontal flight behavior, vertical flight behavior, and acceleration situation; the horizontal flight behavior is specifically the course angle interval between two adjacent moments, which is divided into four intervals at intervals of 90°; the vertical flight behavior is specifically the pitch angle interval between two adjacent moments, including three intervals of descending, climbing, and level flight; the acceleration situation consists of three situations of acceleration, uniform speed, and deceleration between two adjacent moments; if the predicted value of the three-dimensional position coordinates of the non-cooperative UAV B at the t-th moment is only within the conflict protection area of the cooperative UAV A in the horizontal or vertical direction, and the distance between the two in the other direction is less than the flight distance of the UAV running at the maximum speed in a single step, then it is judged whether the probability that the flight behavior interval of the non-cooperative UAV B in the corresponding direction is the same as that of the cooperative UAV A exceeds a predetermined threshold. If it exceeds, it means that there is a conflict risk between the predicted trajectories of the non-cooperative UAV B and the cooperative UAV A at the t-th moment; if it does not exceed, there is no conflict between the cooperative UAV A and the non-cooperative UAV B.

[0143] 23) Based on the results in step 21) and step 22), the prediction label of the conflict behavior between the non-cooperative UAV and the cooperative UAV at the t-th moment is obtained; at the same time, when the relative distances in the horizontal and vertical directions between the true position of the non-cooperative UAV at the t-th moment and the true trajectory coordinates of the cooperative UAV are less than or equal to the set conflict protection area distance, the true label of the conflict event is yes, otherwise it is no.

[0144] 24) Conduct an evaluation on the effectiveness of the conflict identification result.

[0145] 241) Calculate the conflict identification accuracy.

[0146] The conflict identification accuracy is the proportion Acc of the number of correctly predicted samples to the total number of samples. The specific calculation formula is as follows:

[0147]

[0148] In the formula, y l , y' l are the true label and the predicted label of sample l respectively, and n is the number of predicted samples.

[0149] 242) Calculate the precision and recall rate of the conflict identification.

[0150] In the conflict identification of the non-cooperative UAV intrusion scenario, the specific calculation formulas for the precision and recall rate are as follows:

[0151]

[0152]

[0153] Wherein, R is the recognition recall rate, and P is the recognition precision rate;

[0154] 243) Calculation of the comprehensive evaluation index for conflict recognition;

[0155] The emphasis degree of the recognition task on the precision rate P and the recall rate R is characterized by the balance factor α. When α > 1, it emphasizes the recall rate R; when α < 1, it emphasizes the precision rate P; P and R are synthesized into an evaluation index F α , and the specific calculation formula is as follows:

[0156]

[0157] Set the threshold of the comprehensive evaluation index to F α * , when F α > F α * , go to step 3); otherwise, return to step 1).

[0158] In the example, set F α * = 0.8. By predicting the future trajectory of the UAV for 5s, 10s, 15s, and 20s respectively, the conflicts of different time lengths in the intrusion scenario are pre-recognized. Taking the real historical trajectory of the non-cooperative UAV as the training sample, the recognition effect of the UAV conflict recognition is as follows:

[0159] Table 1

[0160] Prediction duration Acc P R <![CDATA[F 0.5 > <![CDATA[F 1.5 > 5s 0.9852 0.9940 0.9880 0.9928 0.9898 10s 0.9396 0.9624 0.9500 0.9598 0.9538 15s 0.9259 0.8972 0.9249 0.9026 0.9162 20s 0.9053 0.8874 0.9114 0.9041 0.9207

[0161] 3) Based on the conflict recognition results in step 2), simulate and simulate the point pair combinations with conflicts to generate Monte Carlo simulation samples;

[0162] Among them, the specific process of step 3) is as follows:

[0163] 31) Establish the kinematic equation of the UAV;

[0164] Establish the kinematic equation of the UAV to ensure the physical feasibility and authenticity of the randomly generated motion trajectory within each adjacent moment. The specific formula is as follows:

[0165]

[0166] Wherein, a is the acceleration vector, V is the velocity vector, T is the thrust vector, K is the empirical drag coefficient, m is the mass of the UAV, and g is the acceleration; among them, K is related to the change of the UAV attitude, and at different tilt angles Under this condition, the relationship between K and the axial resistance coefficient K a , the radial resistance coefficient K s is as follows:

[0167]

[0168] Assume that the drag term of the UAV is axisymmetric, and solve by combining the relevant performance parameters of the UAV in the normal gear and the sports gear to determine the initial thrust and drag. The specific formula is as follows:

[0169]

[0170] In the formula, and are the maximum tilt angles of the UAV in the sports gear and the normal gear respectively, and V1 and V2 are the maximum horizontal speeds of the UAV in the sports gear and the normal gear respectively;

[0171] 32) The thrust vector and attitude inclination angle required when generating a specified velocity vector;

[0172] Assume that the attitude and speed changes occur at the initial stage of each time step. The thrust vector of the UAV is in the same direction as the attitude vector, and the direction of the drag vector is opposite to the direction of the velocity vector. The unit attitude vector and the unit velocity vector The specific calculation formula is as follows:

[0173]

[0174]

[0175] In the formula, ψ and θ are the heading angle and pitch angle respectively, and are the rotation angles around the x-axis and the y-axis (the direction of the UAV's travel) respectively;

[0176] The randomly generated UAV speed and attitude meet the physical feasibility requirements of the motion trajectory. The specific formula is as follows:

[0177]

[0178] In the formula, D is the drag vector. According to the constraint condition When the acceleration vector a = 0, combine Equation (9) and Equation (14), and by iteratively changing and until the thrust vector T is perpendicular to the unit attitude vector Generate the thrust vector and attitude inclination angle required when generating a given velocity vector, so that the randomly generated speed and attitude meet the physical feasibility requirements of the motion trajectory;

[0179] In the example, the performance parameters of the drone are set as shown in Table 2 below:

[0180] Table 2

[0181] Parameter Value UAV mass 242g Maximum horizontal speed (Sport mode - S mode) 16 m / s Maximum tilt angle (Sport mode - S mode) 40° Maximum horizontal speed (Normal mode - P mode) 10 m / s Maximum tilt angle (Normal mode - P mode) 25° Maximum ascent speed 5 m / s Maximum descent speed 3.5 m / s

[0182] 33) Generate Monte Carlo simulation samples;

[0183] For each Monte Carlo simulation sample, the flight trajectory randomly generates the speed, heading angle, and elevation angle at each time step. The initial generation of the speed, heading angle, and elevation angle follows a normal distribution, and Monte Carlo simulation samples of the drone at different times are generated respectively. The specific generation method is as follows:

[0184]

[0185] In the formula, |V| represents the magnitude of the speed; |V target | represents the magnitude of the average speed; |V max | represents the magnitude of the maximum speed; ψ and θ are the heading angle and pitch angle respectively; V target , ψ target , θ target are the means of the speed, heading angle, and pitch angle respectively, and σ V , σ ψ , σ θ are the variances of the speed, heading angle, and pitch angle respectively.

[0186] In the example, set |V max | = 20 m / s, ψ max = 180°, θ max = 90°. Take the target heading angle and pitch angle from the initial position to the prediction moment as the average value of the normal distribution. The target heading angle ψ target and pitch angle θ target respectively point to the traveling direction of the position point at the prediction moment from the current time step. Set the variances of the speed, heading angle, and pitch angle to 0.04, and generate 1000 Monte Carlo simulation sample trajectories at prediction step lengths of 5 seconds, 10 seconds, 15 seconds, and 20 seconds respectively.

[0187] 4) Based on the simulation samples generated in step 3), evaluate the likelihood and severity of the occurrence of dangerous behaviors of non - cooperative drones invading cooperative drones, determine risk assessment indicators, and conduct ratings;

[0188] Among them, the specific process of step 4) is as follows:

[0189] 41) Evaluate the likelihood of the occurrence of dangerous behaviors of non - cooperative drones invading cooperative drones;

[0190] 411) Calculate the predicted invasion probability of the drone;

[0191] The predicted intrusion probability represents the probability of the behavior that a non - cooperative UAV intrudes into the conflict protection area and collision protection area of a cooperative UAV at the prediction moment, and is obtained by calculating the percentage of Monte Carlo simulation samples located in the conflict protection area or collision protection area of the cooperative UAV. The specific calculation formula is as follows:

[0192]

[0193] In the formula, μ1 is the predicted intrusion probability, f is the number of Monte Carlo simulation samples, s k ∈I indicates that the position of the k - th Monte Carlo simulation sample generated is located in the protection area I of the cooperative UAV c where, when c = 1, it is the conflict protection area, and when c = 2, it is the collision protection area;

[0194] In the example, the specific ratings of the predicted intrusion probability are shown in Table 3 below:

[0195] Table 3

[0196] PI rating Parameter value Meaning 1 <![CDATA[0 < μ1 < 0.2]]> Minimum intrusion probability, lowest risk 2 <![CDATA[0.2≤μ1<0.4]]> 3 <![CDATA[0.4≤μ1<0.6]]> 4 <![CDATA[0.6≤μ1<0.8]]> 5 <![CDATA[0.8 ≤ μ1 ≤ 1]]> Maximum intrusion probability, highest risk

[0197] 412) Calculate the predicted conflict time of the UAV;

[0198] When the prediction step is set to i, the predicted position of the cooperative UAV is (x i ,y i ,z i ), and the distribution position of the Monte Carlo simulation samples of the non - cooperative UAV is where, j is the number of generated simulation samples; by calculating the distances from the cooperative UAV to each sample point of the non - cooperative UAV, a sample distance set is obtained, where, is the distance between the predicted point of the cooperative UAV and the p - th simulation sample of the non - cooperative UAV at the prediction step i, and the specific calculation formula is as follows:

[0199]

[0200] The minimum distance between the cooperative UAV and the non - cooperative UAV at the prediction step i The specific calculation formula is as follows:

[0201]

[0202] Assuming that the cooperative UAV is relatively stationary and does not perform any maneuvers, the calculation of the predicted conflict time is carried out under the worst - case intention assumption that the non - cooperative UAV flies straight towards the cooperative UAV at the maximum flight speed. The specific calculation formula is as follows:

[0203]

[0204] where μ2 is the predicted conflict time, is the minimum distance between the cooperative UAV and the non - cooperative UAV at the prediction step i, is the maximum flight speed of the non - cooperative UAV;

[0205] In the example, the specific ratings of the predicted conflict time are shown in Table 4 below:

[0206] Table 4

[0207] TC rating Parameter value(s) Meaning 1 <![CDATA[μ2>20]]> Longest expected time, highest likelihood of collision avoidance, lowest risk 2 <![CDATA[15≤μ2<20]]> 3 <![CDATA[10≤μ2<15]]> 4 <![CDATA[5≤μ2<10]]> 5 <![CDATA[μ2 < 5]]> Shortest expected time, lowest likelihood of collision avoidance, highest risk

[0208] 413) Calculate the predicted conflict intention of the UAV;

[0209] Set the horizontal intention interval of the non - cooperative UAV during the prediction step i → i + 1 as I h , the vertical intention interval I v , and the speed change situation I a When accelerating, I a = 1, when moving at a constant speed, I a = 0, when decelerating, I a = - 1; The predicted position of the cooperative UAV at step i is (x i , y i , z i ); Determine the conflict intention level by combining the relative position of the intention interval and the speed change situation;

[0210] When (x i , y i ) ∈ I h , z i ∈ I v and I a = 1, the non - cooperative UAV will enter the position interval of the cooperative UAV in the horizontal and vertical directions at the next moment, and when accelerating, it means that the non - cooperative UAV has an intrusion intention, and the risk level of an accident is the highest; When and I a = - 1, the position intervals of the non - cooperative UAV in the horizontal and vertical directions are isolated from the cooperative UAV at the next moment, and when decelerating, it means that the non - cooperative UAV maintains a safe separation between UAVs, and the risk level of an accident is the lowest;

[0211] In the example, the specific ratings of the predicted conflict intention are shown in Table 5 below:

[0212] Table 5

[0213]

[0214]

[0215] 414) Calculate the UAV distance avoidance index;

[0216] Use the positions and movement directions between conflicting entities to calculate the UAV distance avoidance index. The specific calculation formula is as follows:

[0217]

[0218] In the formula, μ4 is the distance avoidance index, X is the distance between the cooperative UAV and the non - cooperative UAV at the predicted moment, and D mar is the minimum distance required for the cooperative UAV to safely execute an avoidance maneuver, representing the minimum avoidance range;

[0219] In the example, the specific ratings of the distance avoidance index are shown in Table 6 below:

[0220] Table 6

[0221] DE rating Parameter value Meaning 1 <![CDATA[2 < μ4]]> Longest intrusion distance, highest avoidance safety, lowest risk 2 <![CDATA[1.6 ≤ μ4 < 2]]> 3 <![CDATA[1.3 ≤ μ4 < 1.6]]> 4 <![CDATA[1≤μ4<1.3]]> 5 <![CDATA[μ4≤1]]> Shortest intrusion distance, lowest avoidance safety, highest risk

[0222] 415) Calculate the UAV azimuth avoidance index;

[0223] Use the geometric definition of the velocity obstacle method to determine whether the cooperative UAV and the non - cooperative UAV are on a potentially colliding route. Take the position of the non - cooperative UAV at the predicted moment as the center of a circle, and use the horizontal safety interval of the conflict protection area as the radius to draw a circle; draw two tangent line segments from the position of the cooperative UAV at the predicted moment to form a conical area, which is the collision cone. If the vector sum V h of the velocity V of the cooperative UAV and the velocity V r of the non - cooperative UAV is located within the collision cone, there is a possibility of conflict, and calculate the azimuth avoidance index μ5 as follows: s In the formula, ω is half of the detection angle of the cooperative UAV, V

[0224]

[0225] s represents the vector sum of the velocities of the cooperative UAV and the non - cooperative UAV, and α1 and α2 are the angles formed by the two boundaries of V s and the detection angle of the host;

[0226] In the example, the specific ratings of the azimuth avoidance index are shown in Table 7 below:

[0227] Table 7

[0228]

[0229]

[0230] 416) Determine the probabilities of non - cooperative drones invading the conflict area and the collision area, the expected conflict time between non - cooperative drones and cooperative drones, the expected conflict intention of non - cooperative drones, the distance avoidance index of cooperative drones, the horizontal azimuth avoidance index of cooperative drones, and the vertical azimuth avoidance index of cooperative drones as evaluation indicators for the likelihood of dangerous behavior, and divide them into five levels from 1 to 5 respectively;

[0231] 42) Evaluate the severity of the occurrence of dangerous behavior of non - cooperative drones invading cooperative drones;

[0232] 421) Obtain the drone crash area;

[0233] The drone crash area is divided into a vertical collision area and a horizontal collision area. The specific calculation formula for the vertical collision area is as follows:

[0234] Re v =π(r p +r uav ) 2 (22)

[0235] In the formula, Re v is the vertical collision area, r p is the average radius of the human body, and r uav is the maximum size radius of the drone;

[0236] The calculation formula for the horizontal collision area is as follows:

[0237] Re c =2(r p +r uav )·d+π(r p +r uav ) 2 (23)

[0238] In the formula, Re c is the horizontal collision area, d represents the horizontal distance traveled by the falling drone when it descends to the height of the human body, and d = H p / tanγ; H p is the average height of the human body, and γ is the angle formed by the velocity vector and the horizontal plane or the human body surface impact;

[0239] 422) Estimate the ground collision kinetic energy of the drone;

[0240] Calculate the velocity V imp of the drone collision point, which is obtained through the maximum flight speed V x and the free - fall speed V y relative to the maximum flight height, as follows:

[0241]

[0242] Through V x and V y , solve the collision angle γ of the drone, specifically as follows:

[0243]

[0244] According to the speed magnitude V of the drone at the moment of impact imp and the maximum takeoff mass M, the kinetic energy E of the crashed drone at the impact point can be solved c , specifically as follows:

[0245]

[0246] 423) Calculate the probability of death in a ground collision of the drone;

[0247] Consider the probability of death P under the shielding factor f The specific calculation formula is as follows:

[0248]

[0249] In the formula, P s represents the shielding factor, E h represents the impact energy corresponding to a 50% probability of death when the shielding factor is 6, E d represents the lowest impact energy threshold for causing death, E c is the kinetic energy of the crashed drone at the impact point, μ represents the correction factor, which is used to improve the estimation of low kinetic energy values less than or near the E d threshold, specifically as follows:

[0250]

[0251] In the example, according to the probability of death P f conduct a warning rating of the risk severity, and the specific rating is shown in Table 8 below:

[0252] Table 8

[0253] <![CDATA[P f Rating]]> Parameter value Meaning 1 <![CDATA[P f <0.2]]> Minimum death probability, lowest risk of casualty severity 2 <![CDATA[0.2 ≤ P f <0.4]]> 3 <![CDATA[0.4 ≤ P f <0.6]]> 4 <![CDATA[0.6 ≤ P f < 0.8]]> 5 <![CDATA[P f ≥0.8]]> Maximum death probability, highest risk of casualty severity

[0254] 424) Calculate the number of ground casualties per flight hour;

[0255] Based on formulas (22) - (28), the calculation formula for the number of ground casualties per flight hour of the drone is obtained as follows:

[0256] N p = P b ·Re c ·D p ·P f (29)

[0257] In the formula, N P is the number of ground casualties per flight hour, and P b is the failure probability of the UAV. Since it is considered that conflict or collision events in the intrusion scenario may cause accidents and failures of the UAV, therefore, P b is set to μ1, and Re c is the area of the horizontal collision region, D p is the population density of the UAV operation area, and P f is the probability of death;

[0258] 425) Determine the probability of death and the number of ground casualties per flight hour considering the shielding factor as the evaluation index of the severity of the occurrence of dangerous behaviors, and divide the probability of death into five levels from 1 to 5.

[0259] 5) Based on the risk assessment index determined in step 4), conduct conflict risk assessment and level estimation for the scenario of non-cooperative UAV intrusion into cooperative UAVs;

[0260] Among them, the specific process of step 5) is as follows:

[0261] 51) Set the UAV performance parameters and environmental parameters, including the maximum takeoff mass of the UAV, the maximum operating speed, the maximum size radius, the detection distance, the detection angle, the population density of the commercial area, the shielding coefficient of the commercial area, the average height of the standing human body, the average radius of the human body, and the acceleration due to gravity;

[0262] In the example, the parameter settings of the risk assessment model are specifically as shown in Table 9 below:

[0263] Table 9

[0264] Parameter Parameter value Parameter Parameter value Population density in commercial area <![CDATA[2.41×10 -2 persons / m 2 > Maximum takeoff mass of UAV 3.4 kg Shading coefficient in commercial area 40 Maximum operating speed of UAV 94 km / h Average height of standing human body 1.8m Maximum size radius of UAV 0.3m Average radius of human body 0.25m Detection distance of UAV 0.7～30m Gravitational acceleration <![CDATA[9.8m / s 2 > Detection angle of UAV Horizontal 57°, vertical 44°

[0265] 52) Conduct risk estimation and early warning level determination for the conflict and collision prediction simulation example cases at prediction steps of 5s, 10s, 15s, and 20s respectively, including the probability of intrusion into the conflict and collision protection area, the predicted conflict time, the predicted conflict intention, the distance avoidance index, the azimuth avoidance index in the horizontal and vertical directions, the probability of death, and the number of ground casualties per flight hour.

[0266] In the example, the number of ground casualties per flight hour at prediction steps of 5s, 10s, 15s, and 20s has reached the safe level of 10 -8 , and the other index parameter values and risk level results are as Figure 2 shown.

[0267] The specific application ways of the present invention are numerous. The above are only the preferred embodiments of the present invention. It should be noted 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 assessing the conflict risk of non - cooperative UAVs invading cooperative UAVs, characterized in that, The steps are as follows: 1) Collect real trajectory data of multiple types of drones, process and segment the trajectory data to generate trajectory time series slices; 2) Identify conflicts based on the predicted information of non-cooperative drones and cooperative drones at a future moment and the set protected area range; 3) Based on the conflict identification results in step 2), simulate and simulate the conflicting point pair combinations to generate Monte Carlo simulation samples; 4) Based on the simulation samples generated in step 3), evaluate the possibility and severity of the occurrence of dangerous behaviors of non-cooperative drones invading cooperative drones, determine risk assessment indicators and conduct ratings; 5) Based on the risk assessment indicators determined in step 4), conduct conflict risk assessment and level calculation for the scenario of non-cooperative drones invading cooperative drones; The specific process of step 2) is as follows: 21) Judge the relative position relationship between the predicted position of the non-cooperative drone at a future moment and the conflict protection area of the cooperative drone; 211) Define the conflict protection area and collision protection area of the cooperative drone; Taking the predicted position of the cooperative drone at a future moment as the midpoint of the horizontal and vertical directions of the protection area, considering the two processes of conflict and collision, set the shapes of the conflict protection area and the collision protection area as cylinders, and the conflict protection area is formed by expanding the collision protection area; The size calculation formula of the collision protection area is as follows: Where h is the height of the cylinder in the collision protection area, r is the radius of the cylinder in the collision protection area, h u is the height of the external dimension of the UAV, r u is half of the maximum value of the horizontal dimension of the UAV, G y and G h are the hovering accuracies of the UAV in the vertical and horizontal directions in GPS mode respectively; the cylinder height D ver and radius D hor of the conflict protection area are enlarged by a certain proportion based on the cylinder height and radius of the collision protection area; 212) Judge the conflict interval of the drones; Judge whether the predicted values of the three-dimensional position coordinates of drone B as a particle at the t-th moment are both within the conflict protection area of cooperative drone A in the horizontal and vertical directions. If so, then the following formulas (3) and (4) hold simultaneously, and there is a conflict between drone A and drone B. The specific judgment formulas are as follows: |P(z) i-coo -P(z) i-non | < D ver (4) wherein, P(x) i-coo , P(y) i-coo , P(z) i-coo are respectively the x coordinate in the horizontal direction, the y coordinate in the horizontal direction and the z coordinate in the vertical direction of the cooperative UAV, and P(x) i-non , P(y) i-non , P(z) i-non are respectively the x coordinate in the horizontal direction, the y coordinate in the horizontal direction and the z coordinate in the vertical direction of the non - cooperative UAV; If neither formula (3) nor formula (4) holds, then there is no conflict between drone A and drone B; If only one of formulas (3) and (4) holds, then go to step 22); 22) Based on the judgment result in step 21), if formulas (3) and (4) are not satisfied simultaneously, conflict identification is further carried out in combination with the predicted flight intention of the non-cooperative UAV at future moments; among them, the flight intention within the time period from a certain future moment to the next moment is characterized by three aspects, including horizontal flight behavior, vertical flight behavior, and acceleration situation; the horizontal flight behavior is specifically the course angle interval between two adjacent moments, which is divided into four intervals at intervals of 90°; the vertical flight behavior is specifically the pitch angle interval between two adjacent moments, including three intervals of descent, climb, and level flight; the acceleration situation consists of three situations of acceleration, uniform speed, and deceleration between two adjacent moments; if the predicted value of the three-dimensional position coordinates of the non-cooperative UAV B at the t-th moment is only within the conflict protection area of the cooperative UAV A in the horizontal or vertical direction, and the distance between the two in the other direction is less than the flight distance traveled at the maximum speed under a single UAV step, then it is judged whether the probability that the flight behavior interval of the non-cooperative UAV B in the corresponding direction is in the same interval as the cooperative UAV A exceeds a predetermined threshold. If it exceeds, it means that there is a conflict risk between the predicted trajectories of the non-cooperative UAV B and the cooperative UAV A at the t-th moment; if it does not exceed, there is no conflict between the cooperative UAV A and the non-cooperative UAV B. 23) Based on the results in steps 21) and 22), obtain the prediction label of the conflict behavior between the non-cooperative UAV and the cooperative UAV at the t-th moment; at the same time, when the relative distances in the horizontal and vertical directions between the real position of the non-cooperative UAV at the t-th moment and the real trajectory coordinates of the cooperative UAV are less than or equal to the set conflict protection area distance, the real label of the conflict event is yes, otherwise it is no. 24) Conduct an evaluation on the effectiveness of the conflict identification result. 241) Calculate the conflict identification accuracy. The conflict identification accuracy is the proportion Acc of the number of correctly predicted samples to the total number of samples. The specific calculation formula is as follows: where y l and y l ' are the true label and the predicted label of sample l respectively, and n is the number of predicted samples; 242) Calculate the precision and recall rate of the conflict identification. In the conflict identification of the non-cooperative UAV intrusion scenario, the specific calculation formulas for the precision and recall rate are as follows: In the formula, R is the recall rate of identification, and P is the precision rate of identification. 243) Calculate the comprehensive evaluation index of the conflict identification. The emphasis degree of the recognition task on the precision P and the recall R is characterized by the balance factor α. When α > 1, it emphasizes the recall R; when α < 1, it emphasizes the precision P. P and R are combined into an evaluation index F α , and the specific calculation formula is as follows: Set the threshold of the comprehensive evaluation index to F α * , when F α >F α * , go to step 3); otherwise, return to step 1).

2. The conflict risk assessment method for non - cooperative drones invading cooperative drones according to claim 1, wherein, The specific process of step 1) is as follows: Collect real trajectory data of various types of UAVs. Each trajectory point in the trajectory data includes timestamp, identity ID, latitude, longitude, relative flight altitude, speed, and attitude angle information. Convert the longitude, latitude, and altitude in the WGS-84 coordinate system into (x, y, z) in the Cartesian coordinate system. Perform equidistant sampling on the real trajectory data of the UAVs, use the sliding time window method to slice and divide the trajectory sequence data of the UAVs, and perform normalization processing to generate the trajectory time series slices of various types of UAVs.

3. The conflict risk assessment method for non-cooperative drones invading cooperative drones according to claim 1, wherein The specific process of step 3) is as follows: 31) Establish the UAV kinematic equation. Establish the UAV kinematic equation to ensure the physical feasibility and authenticity of the randomly generated motion trajectories within each adjacent moment. The specific formula is as follows: Wherein, a is the acceleration vector, V is the velocity vector, T is the thrust vector, K is the empirical drag coefficient, m is the mass of the UAV, and g is the acceleration; among them, K is related to the change of the UAV attitude, and at different inclination angles , the relationship between K and the axial drag coefficient K a , the radial drag coefficient K s is: Set the drag term of the UAV to be axisymmetric, and solve it by combining the relevant performance parameters of the UAV in the normal gear and the sports gear to determine the initial thrust and drag. The specific formula is as follows: Wherein, and are respectively the maximum tilt angles of the drone in the motion gear and the normal gear, and V1 and V2 are respectively the maximum horizontal speeds of the drone in the motion gear and the normal gear; 32) Generate the thrust vector and attitude inclination angle required when generating a specified velocity vector; Set the attitude and speed changes to occur at the beginning of each time step. The thrust vector of the UAV is in the same direction as the attitude vector, and the direction of the drag vector is opposite to the direction of the speed vector. The unit attitude vector and the unit speed vector The specific calculation formula is as follows: where ψ and θ are the heading angle and pitch angle respectively, and are the rotation angles about the x-axis and y-axis respectively; The randomly generated UAV speed and attitude meet the physical feasibility requirements of the flight trajectory. The specific formula is as follows: where D is the drag vector, according to the constraint conditions In the case where the acceleration vector a = 0, combining Equation (9) and Equation (14), by iteratively changing and until the thrust vector T and the unit attitude vector are perpendicular, the thrust vector and the attitude inclination angle required to generate a given velocity vector are generated, so that the randomly generated velocity and attitude meet the physical feasibility requirements of the motion trajectory; 33) Generate Monte Carlo simulation samples; For each Monte Carlo simulation sample, the flight trajectory randomly generates speed, heading angle, and elevation angle at each time step. The initial generation of speed, heading angle, and elevation angle follows a normal distribution, and Monte Carlo simulation samples of the UAV at different times are generated respectively. The specific generation method is as follows: where |V| represents the magnitude of the velocity; |V target | represents the magnitude of the average velocity; |V max | represents the magnitude of the maximum velocity; ψ and θ are the course angle and the pitch angle respectively; V target 、ψ target 、θ target are the means of the velocity, the course angle and the pitch angle respectively, and σ V 、σ ψ 、σ θ are the variances of the velocity, the course angle and the pitch angle respectively.

4. The conflict risk assessment method for non - cooperative drones invading cooperative drones according to claim 1, wherein The specific process of step 4) is as follows: 41) Evaluate the possibility of a non - cooperative UAV invading a cooperative UAV's dangerous behavior; 411) Calculate the predicted invasion probability of the UAV; The predicted invasion probability represents the probability of the non - cooperative UAV invading the conflict protection area and collision protection area of the cooperative UAV at the prediction moment, and is obtained by calculating the percentage of Monte Carlo simulation samples located in the conflict protection area or collision protection area of the cooperative UAV. The specific calculation formula is as follows: where μ1 is the predicted intrusion probability, f is the number of Monte Carlo simulation samples, and s k ∈I indicates that the position of the k-th generated Monte Carlo simulation sample is within the protected area I of the cooperative UAV, where c when c = 1, it is a conflict protected area, and when c = 2, it is a collision protected area; 412) Calculate the predicted conflict time of the UAV; When the prediction step is set to i, the predicted position of the cooperative UAV is (x i , y i , z i ). The Monte Carlo simulation sample distribution position of the non-cooperative UAV is where j is the number of generated simulation samples; by calculating the distances from the cooperative UAV to each sample point of the non-cooperative UAV, the sample distance set is the distance between the predicted point of the cooperative UAV and the p-th simulation sample of the non-cooperative UAV at the prediction step i. The specific calculation formula is as follows: The minimum distance between the cooperative UAV and the non - cooperative UAV at prediction step i The specific calculation formula is as follows: Assume that the cooperative UAV is relatively stationary and does not perform any maneuvers, and calculate the predicted conflict time based on the worst - case scenario that the non - cooperative UAV flies straight towards the cooperative UAV at the maximum flight speed. The specific calculation formula is as follows: where μ2 is the predicted conflict time, is the minimum distance between the cooperative UAV and the non - cooperative UAV at the prediction step i, is the maximum flight speed of the non - cooperative UAV; 413) Calculate the predicted conflict intention of the UAV; Set the horizontal intention interval of the non - cooperative UAV in the prediction step from \(i\) to \(i + 1\) as \(I\). h , the vertical intention interval \(I\). v , the speed change situation \(I\). a When accelerating, \(I\). a = 1, when moving at a constant speed, \(I\). a = 0, when decelerating, \(I\). a = - 1; The predicted position of the cooperative UAV at step \(i\) is \((x\). i , y\). i , z\). i ); Judge the conflict intention level by combining the relative position of the intention interval and the speed change situation. When (x i , y i ) ∈ I h , z i ∈ I v and I a = 1, the non - cooperative UAV will enter the position interval of the cooperative UAV in the horizontal and vertical directions at the next moment, and when accelerating, it indicates that the non - cooperative UAV has an intrusion intention and the risk level of an accident is the highest; when and I a = - 1, the position intervals of the non - cooperative UAV in the horizontal and vertical directions are isolated from those of the cooperative UAV at the next moment, and when decelerating, it indicates that the non - cooperative UAV maintains a safe interval between UAVs and the risk level of an accident is the lowest; 414) Calculate the distance avoidance index of the UAV; Use the positions and motion directions between the conflict entities to calculate the distance avoidance index of the UAV. The specific calculation formula is as follows: In the formula, μ4 is the distance avoidance index, X is the distance between the cooperative UAV and the non - cooperative UAV at the predicted moment, and D mar is the minimum distance required for the cooperative UAV to safely execute the avoidance maneuver, representing the minimum avoidance range; 415) Calculate the azimuth avoidance index of the UAV; Use the geometric definition of the speed obstacle method to determine whether the cooperative UAV and the non-cooperative UAV are on a potential collision course. Take the position of the non-cooperative UAV at the prediction moment as the center of the circle, and the horizontal safety interval of the conflict protection area as the radius to draw a circle; take two tangent line segments of the circle with the position of the cooperative UAV at the prediction moment, and the conical area formed is the collision cone; if the vector sum V h of the speed V r of the cooperative UAV and the speed V s of the non-cooperative UAV is located within the collision cone, there is a possibility of conflict, and calculate the azimuth avoidance index μ5 as follows: where ω is half of the detection angle of the cooperative UAV, and V s represents the sum of the velocity vectors of the cooperative UAV and the non - cooperative UAV, and α1 and α2 are the angles formed by the two boundaries of V s with the detection angle of the host respectively; 416) Determine the probability of the non - cooperative UAV invading the conflict area, the probability of the non - cooperative UAV invading the collision area, the predicted conflict time between the non - cooperative UAV and the cooperative UAV, the predicted conflict intention of the non - cooperative UAV, the distance avoidance index of the cooperative UAV, the horizontal azimuth avoidance index of the cooperative UAV, and the vertical azimuth avoidance index of the cooperative UAV as the evaluation indicators for the possibility of dangerous behavior occurring; 42) Evaluate the severity of the non - cooperative UAV invading the cooperative UAV's dangerous behavior; 421) Obtain the UAV crash area; The UAV crash area is divided into a vertical collision area and a horizontal collision area. The specific calculation formula for the vertical collision area is as follows: Re v = π(r p + r uav ) 2 (22) where Re v is the vertical collision area, r p is the average radius of the human body, and r uav is the maximum size radius of the drone; The calculation formula for the horizontal collision area is as follows: Re c = 2(r p + r uav )·d + π(r p + r uav ) 2 (23) Where Re c is the horizontal collision area, d represents the horizontal distance traveled by the falling drone when it descends to the height of the human body, and d = H p / tanγ; H p is the average height of the human body, and γ is the angle formed by the velocity vector and the horizontal plane or the human body surface impact; 422) Estimate the ground collision kinetic energy of the UAV; Calculate the collision point velocity V of the drone imp , through the maximum flight speed V x and the free fall velocity V relative to the maximum flight height y obtained as follows: Through V x and V y , solve the collision angle γ of the drone as follows: According to the speed V of the drone at the moment of impact imp and the maximum take-off mass M, the kinetic energy E of the crashed drone at the impact point can be solved c as follows: 423) Calculate the ground collision fatality probability of the UAV; Probability of death \(P\) considering the shielding factor f The specific calculation formula is as follows: where P s represents the shielding factor, E h represents the impact energy corresponding to a 50% probability of death when the shielding factor is 6, E d represents the minimum impact energy threshold for causing death, E c is the kinetic energy of the crashed UAV at the impact point, and μ represents a correction factor used to improve the estimation of low kinetic energy values less than or near the E d threshold, as follows: 424) Calculate the number of ground casualties per flight hour; Based on formula (22) - formula (28), the calculation formula for the number of ground casualties per flight hour of the UAV is as follows: N p = P b · Re c · D p · P f (29) Where N P is the number of ground casualties per flying hour, P b is the UAV failure probability. Since conflict or collision events in the intrusion scenario may cause accidents and failures of the UAV, therefore, P b is set to μ1, Re c is the horizontal collision area, D p is the population density of the UAV operation area, and P f is the probability of death; 425) Determine the fatality probability considering the shielding factor and the number of ground casualties per flight hour as the evaluation indicators for the severity of dangerous behavior occurring.

5. The conflict risk assessment method for a non - cooperative drone invading a cooperative drone according to claim 1, wherein, The specific process of step 5) is as follows: 51) Set the performance parameters and environmental parameters of the UAV, including the maximum takeoff mass of the UAV, the maximum operating speed, the maximum size radius, the detection distance, the detection angle, the population density of the commercial area, the shielding coefficient of the commercial area, the average height of standing humans, the average radius of the human body, and the acceleration due to gravity; 52) Calculate the risk and determine the warning level for the conflict and collision prediction simulation example cases at prediction time steps of 5 s, 10 s, 15 s, and 20 s respectively, including the probability of intrusion into the conflict and collision protection area, the expected conflict time, the expected conflict intention, the distance avoidance index, the azimuth avoidance index in the horizontal and vertical directions, the probability of death, and the number of ground casualties per flight hour.

Images