Method for calculating safety interval between unmanned aerial vehicle and manned aircraft

Abstract

The invention discloses a safety interval calculation method for an unmanned aerial vehicle and a manned aircraft, and relates to the technical field of aviation safety. The safe interval calculation method comprises the following steps: step 1, firstly obtaining the length, wingspan and height of a first target, then obtaining the length, wingspan and height of a second target, marking the first target as a mass point at the moment, and making up for an error caused by volume loss when the first target is regarded as the mass point. A new collision model is used for carrying out quantitative analysis on the collision probability of a manned aerial vehicle and an unmanned aerial vehicle so as to find the relation between the flight interval and the collision probability and calibrate a reasonable safety interval, and the requirement of the safety interval with a fixed value for the airspace cost is high. Therefore, the dynamic safety interval is found by adopting interval division iterative calculation, and a more flexible dynamic safety interval between the unmanned aerial vehicle and the manned vehicle is given on the basis, so that the construction of a fusion flight system is accelerated.

Description

Technical Field

[0001] The present invention relates to the technical field of aviation safety, and specifically to a method for calculating the safety interval between an unmanned aerial vehicle and a manned aircraft. Background Art

[0002] Air transportation has been rapidly popularized. Due to the rise of intelligent logistics technology, air logistics carried by fixed-wing unmanned aerial vehicles has started trial operation. Restricted by many conditions such as the number of airports and the scope of civil airspace, the operating contradictions between civil airliners and unmanned aerial vehicles are relatively prominent, especially the contradiction in airspace flight safety. In recent years, many experts and scholars have conducted a large number of studies on the collision risk of aircraft in the airspace and the conflict resolution of aircraft.

[0003] The existing unmanned aerial vehicle approximates the collision box as a sphere, and then calculates the collision probability by approaching a spherical unmanned aerial vehicle with position error to another spherical collision box. The advantage of these regular geometric models is that they reduce unnecessary redundant calculations, but they do not fit well with real aircraft, and do not fully consider the volume and positioning error of the aircraft. Therefore, the present invention proposes a method for calculating the safety interval between an unmanned aerial vehicle and a manned aircraft to solve the above-mentioned problems. Summary of the Invention

[0004] (1) Technical Problems to be Solved

[0005] In view of the deficiencies of the prior art, the present invention provides a method for calculating the safety interval between an unmanned aerial vehicle and a manned aircraft, which solves the problems mentioned in the above background art.

[0006] (2) Technical Solutions

[0007] To achieve the above objectives, the present invention is realized through the following technical solutions: A method for calculating the safety interval between an unmanned aerial vehicle and a manned aircraft includes the following steps:

[0008] Step 1: First, obtain the length, wingspan, and height of the first target, and then obtain the length, wingspan, and height of the second target. At this time, mark the first target as a mass point, and compensate for the error caused by the volume loss of regarding the first target as a mass point. Establish a collision model based on the first target and the second target, where the first target is specifically an unmanned aerial vehicle and the second target is specifically a manned aircraft;

[0009] Step 2: When GPS is used for positioning and navigation, obtain the multipath effect and IMU drift generated during the communication process between the first target and the second target, and then obtain the trajectory of the velocity component changing with time. When obtaining the velocity in each direction, the actual velocity vector can be obtained through the method of vector synthesis.

[0010] Step 3: Add a three-dimensional Gaussian distribution adjoint error in the direction of the first target velocity vector and a Laplace distribution adjoint error in the direction of the second target velocity vector. Then, consider the first target as a particle approaching the collision box in three-dimensional space, and calculate the volume of the first target within the collision box through multiple integrals to represent the collision probability, thereby obtaining the relationship between the interval and the collision probability.

[0011] The specific method for compensating the error caused by the volume loss of considering the first target as a particle in Step 1 is as follows:

[0012] Establish a collision box in combination with the size of the UAV fuselage.

[0013] S i <<λ i i∈x,z j∈x,y (1)

[0014] Based on the classical spherical and circular cuboid collision boxes, the collision box parameters with the UAV volume are defined as follows in this paper:

[0015]

[0016] Then, in the three-dimensional coordinate system, the ellipsoidal collision box is determined by formula (5):

[0017]

[0018] Taking the center of gravity of the airliner as the centroid and the ellipsoid with parameters abh as the collision box, when the UAV with a particle flying at a certain velocity vector approaches the large aircraft with a collision box flying at a certain velocity vector in the integrated airspace, it is regarded as the occurrence of a collision risk. When the UAV moves in the joint airspace for time T, its movement trajectory is ΩT, and the collision box area attached to the large aircraft is ΦT. The size of the collision probability is determined by the size of the intersection of ΩT and ΦT:

[0019] C=Φ T ∩Ω T (6).

[0020] Preferably, the specific method for then obtaining the trajectory of the velocity component changing with time and obtaining the actual velocity vector through vector synthesis when obtaining the velocities in each direction in Step 2 is as follows:

[0021] In the positioning principle, it is known that when the GPS signal is affected by multiple effects such as superposition shadows, its positioning error is transformed into a Gaussian process. The error ε generated by the UAV follows a Gaussian normal distribution with a mean of u,x = 0 and a standard deviation of δ=(4*16.12)*0.5 = 0.322 km in the running speed direction [9], that is, formula (7):

[0022]

[0023] For large aircraft, the main positioning error mainly comes from navigation accuracy and meteorological conditions. According to the required performance navigation principle, within the range of n nautical miles during flight, there is a 95% probability that large aircraft deviate from the planned flight path, and its deviation f(∈) follows a Laplace distribution with an expectation of 0

[10] , that is, Equation (8):

[0024]

[0025] To eliminate the influence of positioning error, in this study, a Gaussian distribution is added to the UAV speed vector, and a Laplace distribution is added to the direction of the manned aircraft speed vector for simulation calculation. When the speed vector of an object is v, and the position along the speed direction follows

[0026] When starting from the Laplace, a stochastic differential equation can be used to describe the change of the speed component. Assume that the speed vector of the object is: v = (vx, vy, vz), where vx, vy, and vz respectively represent the speed components on the xyz axes. Consider the differential equation in Equation (9):

[0027] dv = -αvdt + δdW (9)

[0028] dv represents the small change of the speed component, α is the attenuation coefficient, v is the speed component, dt is the time element, δ is the diffusion coefficient, and dW is the differential of the Wiener process (or Brownian motion). This equation contains two terms: the attenuation term and the diffusion term. The attenuation term (-ɑv) makes the speed component tend to a stable mean value. The diffusion term (δdW) introduces random noise, causing the speed component to fluctuate around its mean value. dW is the differential of the Wiener process, representing random changes. The solution of this stochastic differential equation can be simulated by the Euler numerical method, using the discretized time step Δt to approximate the differential equation to simulate the change of the speed component over time to obtain the trajectory of the speed component.

[0029] Assume that the speed vector v is expressed as v = (vx, vy, vz), and these speed components follow a Gaussian distribution along the speed direction. Then, the change of the speed component can be expressed as Equations (10 - 12):

[0030] dv x = -α x v x dt + δ x dW x (10

[0031] dv y = -α y v y dt + δ y dW y (11

[0032] dv z = -α z v z dt + δ z dW z (12

[0033] dvx, dvy, and dvz represent the small changes in the velocity components along the x, y, and z axes respectively. ɑx, ɑy, and ɑz are the corresponding attenuation coefficients, vx, vy, and vz are the velocity components, δx, δy, and δz are the corresponding diffusion coefficients, and dWx, dWy, and dWz are the differentials of the Wiener process. This formula includes the changes in the velocity components in three dimensions. Each velocity component is affected by an attenuation term and a diffusion term. The attenuation term makes the velocity component tend to a stable mean value, and the diffusion term introduces random noise, causing the velocity component to fluctuate around its mean value.

[0034] Preferably, in the third step, the volume of the first target within the collision box is calculated by double integration to represent the collision probability. The specific way to obtain the relationship between the interval and the collision probability is as follows:

[0035] Suppose the position vectors of the UAV and the manned aircraft in the fusion airspace are rH and rM respectively, and the velocity vectors are Vm and Vh respectively. Since the flight speed of the UAV is relatively small compared to that of large aircraft, for the convenience of subsequent calculations, the process of the UAV approaching the large aircraft is regarded as a process of flying at a constant speed v. Due to the variable speed error of the UAV, this speed cannot be simply set as the theoretical average speed of the UAV. If the stable flight speed is v when the aircraft approaches, according to the principle of speed obstacle method for obstacle avoidance, when the speed of the UAV changes from v to v1 within tf, it satisfies formula (13)

[11]

[0036]

[0037] amax is the maximum acceleration that conforms to the characteristics of the UAV. Affected by technical limitations and the frequency of crosswinds in the fusion airspace, amax satisfies a one-dimensional Gaussian distribution with a mean of amax and a variance of 0.01:

[0038]

[0039] If the maximum speed of the UAV during the overall stable flight stage is Vmax and the minimum speed is Vmin, then from formula (15), the stable flight speed Vau of the UAV considering the error can be obtained:

[0040]

[0041] Within time T, the relative position vector Δr and relative velocity vector Δv between the UAV and the airliner can be expressed as equations (16) and (17) respectively:

[0042] Δr e =r H -r M (16)

[0043] Δv e =v H -v M (17)

[0044]

[0045] The relative nominal position of the UAV and the airliner is equation (18); for probability calculation, the aircraft needs to be transformed from the geodetic coordinates to the body coordinates. According to the coordinate transformation theorem and the Euler angle theorem, the body coordinate system can be obtained by left-multiplying the geodetic coordinate system with the coordinate transformation matrix P

[14]

[0046]

[0047] Through coordinate transformation, the relative position vector re is expressed as formula (20):

[0048] Δr e =r H -r M =P H r H -P M r M (20)

[0049] PM is the transformation matrix from the body coordinate system of the target UAV to the global coordinate system; PH is the transformation matrix from the body coordinate system of the airliner to the global coordinate system. Given the relative position vector and relative velocity vector in the body coordinate system, assuming that the flight trajectories of the UAV and the large aircraft are independent of each other and both follow a one-dimensional normal distribution, then the relative position vector △re also follows a Gaussian distribution

[15] , that is, formula (21):

[0050]

[0051] When the UAV and the aircraft are in relative motion, the collision probability can be regarded as the integral ratio of a probability density function that follows a three-dimensional Gaussian distribution within the collision box area. Let D be the motion area of the UAV within the collision box, then D satisfies formula (22):

[0052]

[0053] Where p is used to represent the three-dimensional space vector within the collision sphere of the airliner

[16] .

[0054] Suppose that at time T1, the point UAV is located at A1 within D. At this time, the relative position vector between the manned aircraft and the UAV is r1. An ellipsoidal surface that takes the center of gravity of the manned aircraft as the center of the sphere, passes through point A1, and is three-dimensionally similar to the ellipsoidal collision box surface is defined as the collision probability surface S1 at this moment. The numerical value of the surface integral of the probability density function over this collision probability surface is defined as PT1, that is, Equation (23):

[0055]

[0056] At time Tn, it is located at An. At this time, the relative position vector between the manned aircraft and the UAV is rn. A spherical surface with the center of gravity of the manned aircraft as the center of the circle and |rn| as the radius is defined as the collision surface Sn at this moment. At this time, PTn is determined by Equation (24):

[0057]

[0058] If the UAV has experienced a total of T1, T2... Tn points within T, then the collision surface P at this time is expressed as Formula (25):

[0059]

[0060] Since the UAV is a continuously changing process within the collision box, T is continuous. Therefore, when the UAV's position changes from A1 to Tn within the collision box during T, its collision volume is transformed into a triple integral, that is, Formula (26):

[0061]

[0062] D1 is the collision surface area region traversed by the UAV during the time T. The further collision probability can be obtained from the following Formula (27):

[0063]

[0064] where DM is the volume of the collision box.

[0065] The integration result of the above formula can be obtained by the method of numerical integration. The above formula gives the calculation method of the collision probability, but does not give the inner collision surface nodes and outer collision surface nodes of the UAV's movement within the collision box. Therefore, a specific numerical value needs to be found within the obtained change interval to represent the maximum collision probability. Suppose the current time is, and the initial position vectors of the UAV and the manned aircraft are ra and rb respectively. The relative position vector can be expressed as Formula (28):

[0066] r0 = r a -r b (28)

[0067] Let the relative velocity between the UAV and the airliner be Vr. Since the velocity can be considered constant at the moment approaching the collision area, the relative position vector of the two aircraft at this time can be expressed by Equation (29):

[0068] r (t) =r0+v r (t) (29)

[0069] In the formula, r0 is the relative position vector of the two aircraft. Let ρ(t) = r2(t). To find its extreme value, its derivative function can be set to 0, that is, Equation (30):

[0070]

[0071] Suppose the solution of the above formula is tume. When t = tume, the position vector must take the minimum value. At this time, the double integral within the region defined by the position of the UAV and the outer boundary of the collision box is the maximum collision probability

[18]

[19] .

[0072]

[0073] The minimum value of the relative position vector can be calculated through Equations (31) and (32), and then the maximum collision probability between the UAV and the manned aircraft at a certain relative position can be calculated through Equation (27).

[0074] (III) Beneficial Effects

[0075] The present invention provides a method for calculating the safety interval between a UAV and a manned aircraft. Compared with the prior art, it has the following beneficial effects:

[0076] The method for calculating the safety interval between the UAV and the manned aircraft quantifies and analyzes the collision probability between the manned aircraft and the UAV by using a new collision model, aiming to find the relationship between the flight interval and the collision probability so as to calibrate a reasonable safety interval. Since a fixed-value safety interval requires a high cost for airspace, this article uses interval division and iterative calculation to find the dynamic safety interval range, and on this basis, gives a more flexible dynamic safety interval between the UAV and the manned aircraft to accelerate the construction of the integrated flight system. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 is a cuboid collision template shown in an embodiment of the present application; Figure 2 is a spherical collision template shown in an embodiment of the present application; Figure 3 is an improved collision template shown in an embodiment of the present application; Figure 4 is a collision schematic diagram shown in an embodiment of the present application; Figure 5Relationship between longitudinal spacing and collision probability shown in the embodiments of this application; Figure 6 Area between the interval of 500m and the longitudinal spacing axis shown in the embodiments of this application; Figure 7 Obvious change interval of longitudinal spacing collision probability shown in the embodiments of this application; Figure 8 Collision probability change when d = 250m shown in the embodiments of this application; Figure 9 Collision probability change when d = 125m shown in the embodiments of this application; Figure 10 Collision probability change when d = 62.5m shown in the embodiments of this application; Figure 11 Collision probability change when d = 31.25m shown in the embodiments of this application; Figure 12 Area between the interval of 100m and the vertical spacing axis shown in the embodiments of this application; Figure 13 Obvious change interval of vertical spacing collision probability shown in the embodiments of this application; Figure 14 Collision probability change when d = 50m shown in the embodiments of this application; Figure 15 Collision probability change when d = 25m shown in the embodiments of this application; Figure 16 Collision probability change when d = 12.5m shown in the embodiments of this application; Figure 17 Flowchart of the safety interval calculation method for unmanned aerial vehicles and manned aircraft shown in the embodiments of this application. Detailed implementation manners

[0078] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0079] Please refer to Figure 1 , the embodiments of the present invention provide a technical solution:

[0080] A safety interval calculation method for unmanned aerial vehicles (UAVs) and manned aircraft. The Reich model is regarded as a classic model in aircraft collision models. It approximates an aircraft as a collision box with a geometric shape, and adds a critical layer for calculating the collision probability around the collision box. The collision probability between aircraft is characterized by the overlapping direction and size of the critical layer in space. The parameters of the collision box are generally a cuboid or a hexahedron similar to a cuboid with λi, i ∈ (x, y, z). λx, λy, and λz are the length, wingspan, and height of the aircraft fuselage respectively. The parameters of the critical layer Sj, j ∈ (x, y, z) are determined by the lateral and longitudinal intervals. [7] 。

[0082] The Event model suitable for calculating the collision probability in modern navigation removes the critical layer. Subsequently, a large number of scholars have optimized it based on the Event model. In 2022, Professor Zhang Honghai approximated the collision box as a sphere according to the characteristics of autogyro UAVs, and then calculated the collision probability by bringing a spherical UAV with position error closer to another spherical collision box. [8] The advantage of these regular geometric models is that they reduce unnecessary redundant calculations, but they do not match the real aircraft well and do not fully consider the volume and positioning error of the aircraft. To solve these problems, this paper proposes a calculation method based on an improved Event model.

[0083] 1.2 Model establishment

[0084] Considering the shape of the aircraft

[0085] Make the following assumptions:

[0086] 1. When not considering the positioning error, the collision template is defined as a regular ellipsoid that conforms to the Event template parameters.

[0087] 2. There is no obvious sudden change in speed during the flight of UAVs and manned aircraft.

[0088] 3. The frequencies of headwinds and tailwinds received by the UAV during flight are equal, there is no flight error caused by wind speed, and there is no extreme weather.

[0089] For the convenience of calculation, this paper regards the UAV as a particle. Suppose the length, wingspan, and height of a large manned aircraft in the joint airspace are λx, λy, and λz respectively, and the length, wingspan, and height of the UAV are Sx, Sy, and Sz respectively. According to the common UAVs and large aircraft operating in the current airspace, their lengths, wingspans, and heights satisfy formula (1) but do not satisfy Sz < λz. Therefore, the volume of the UAV cannot be ignored when the UAV ascends and the large aircraft descends. In the actual flight process, both aircraft have components in the z direction. To compensate for the error caused by the volume loss of regarding the UAV as a particle, a collision box needs to be established in combination with the UAV body size.

[0090] S i <<λ i i∈x,z j∈x,y (1)

[0091] Table 1 UAV parameters

[0092]

[0093] Table 2 Large civil aircraft parameters

[0094]

[0095] According to the classical spherical and circular cuboid collision boxes, the collision box parameters with the UAV volume are defined as follows:

[0096]

[0097] Then, in the three-dimensional coordinate system, the ellipsoidal collision box is determined by formula (5):

[0098]

[0100] Taking the center of gravity of the airliner as the centroid and the ellipsoid with parameters a, b, and h as the collision box, when flying in the fusion airspace, when the UAV given to the particle flies with a certain velocity vector and the large aircraft given to the collision box flies with a certain velocity vector and approaches each other, it is regarded as the occurrence of a collision risk. When flying in the joint airspace for T, the movement trajectory of the UAV is ΩT, and the collision box area attached to the large aircraft is ΦT. The size of the collision probability is determined by the size of the intersection of ΩT and ΦT:

[0101] C=Φ T ∩Ω T (6)

[0103] 2 Positioning error

[0104] The positioning error of the drone is mainly caused by the multipath effect, IMU drift, and signal frequency dispersion during the communication process when GPS is used for positioning and navigation. The multipath effect and satellite clock error are the largest error sources during the flight of the drone. The multipath effect occurs because scattered waves are generated when radio waves encounter various obstacles during transmission, resulting in the distortion of the superimposed signal of the positioning information received at the receiving end. Frequency dispersion occurs in its signal frequency domain, and the fading of the signal envelope is usually a stationary narrowband Gaussian process, that is, Rayleigh fading. In the positioning principle, it can be seen that when GPS signals are affected by multiple effects such as shadow superposition, the positioning error is transformed into a Gaussian process. The error ε of the drone follows a Gaussian normal distribution with a mean of u,x = 0 and a standard deviation of δ=(4*16.12)*0.5 = 0.322 km in the direction of the running speed. [9] That is, Equation (7):

[0105]

[0106] For large aircraft, its main positioning error comes from navigation accuracy and meteorological conditions. According to the required performance navigation principle, within the range of n nautical miles during flight, there is a 95% probability that large aircraft deviate from the planned flight path, and its deviation f(∈) follows a Laplace distribution with an expectation of 0.

[10] That is, Equation (8):

[0107]

[0108] To eliminate the influence of positioning errors, in this study, a Gaussian distribution is added to the drone's velocity vector, and a Laplace distribution is added to the velocity vector direction of the manned aircraft for simulation calculations. When the velocity vector of an object is v and its position along the velocity direction follows

[0109] When starting from the Laplace, the change of the velocity component can be described by a stochastic differential equation. Assume that the velocity vector of the object is: v=(v x ,v y ,v z ), where v x ,v y ,v z represent the velocity components on the xyz axes respectively. Consider the differential equation in Equation (9):

[0110] dv = -αvdt + δdW (9)

[0111] The dv represents a small change in the velocity component, α is the attenuation coefficient, v is the velocity component, dt is the time element, δ is the diffusion coefficient, and dW is the differential of the Wiener process (or Brownian motion). This equation contains two terms: the attenuation term and the diffusion term. The attenuation term (-ɑv) makes the velocity component tend to a stable mean value. The diffusion term (δdW) introduces random noise, causing the velocity component to fluctuate around its mean value. dW is the differential of the Wiener process, representing random variation. The solution of this stochastic differential equation can be simulated by the Euler numerical method, using a discretized time step Δt to approximate the differential equation to simulate the change of the velocity component over time to obtain the trajectory of the velocity component.

[0112] Assume the velocity vector v is expressed as v = (v x , v y , v z ). These velocity components follow a Gaussian distribution along the velocity direction. Then, the change in the velocity component can be expressed by the formulas (10)-(12):

[0113] dv x = -α x v x dt + δ x dW x (10

[0114] dv y = -α y v y dt + δ y dW y (11

[0115] dv z = -α z v z dt + δ z dW z (12

[0116] dv x 、dv y and dv z represent the small changes in the velocity components on the x, y, and z axes respectively. ɑ x 、ɑ y and ɑ z are the corresponding attenuation coefficients, v x 、v y and v z are the velocity components, δx, δy, and δz are the corresponding diffusion coefficients, and dW x 、dW y and dW zIt is the differential of the Wiener process. This formula includes the variations of the velocity components in three dimensions. Each velocity component is affected by a decay term and a diffusion term. The decay term makes the velocity component tend to a stable mean value, and the diffusion term introduces random noise, causing the velocity component to fluctuate around its mean value. By integrating this formula, the trajectory of the velocity component over time can be obtained. When obtaining the velocities in each direction, the actual velocity vector can be obtained through vector synthesis.

[12]

[13] .

[0117] 3 Collision probability calculation

[0118] Since the position errors of the aircraft follow Gaussian distribution and Laplace distribution respectively, in this paper, a three-dimensional Gaussian distribution adjoint error is added in the direction of the UAV velocity vector, and a Laplace distribution adjoint error is added in the direction of the large aircraft velocity vector. Then, the UAV is regarded as a particle approaching the collision box in three-dimensional space, and the volume of the UAV in the collision box is calculated through double integral to represent the collision probability, so that the relationship between the interval and the collision probability can be obtained.

[0119] Suppose the position vectors of the UAV and the manned aircraft in the fusion airspace are r H r M The velocity vectors are V m and V h , because for the flight speed of large aircraft, the flight speed of the UAV is relatively small. Therefore, for the convenience of subsequent calculations, the process of the UAV approaching the large aircraft is regarded as a process of flying at a constant speed v. Due to the variable speed error of the UAV, this speed cannot be simply set as the theoretical average speed of the UAV. If the steady flight speed is v when the aircraft approaches, according to the principle of velocity obstacle method for obstacle avoidance, when the speed of the UAV changes from v to v1 within tf, it satisfies formula (13)

[11]

[0120]

[0121] a max is the maximum acceleration that combines the characteristics of the UAV itself. Affected by technical limitations and the frequencies of crosswinds in the fusion airspace, a max satisfies a one-dimensional Gaussian distribution with a mean of a max and a variance of 0.01:

[0122]

[0123] If the maximum speed of the UAV during variable speed in the overall steady flight stage is V max and the minimum speed is V min , then from formula (15), the steady flight speed V au of the UAV considering errors can be obtained:

[0124]

[0125] Within time T, the relative position vector Δr and relative velocity vector Δv between the UAV and the airliner can be respectively expressed as in Equations (16) and (17):

[0126] Δr e = r H - r M (16)

[0127] Δv e = v H - v M (17)

[0128]

[0129] The relative nominal position of the UAV and the airliner is in Equation (18); for probability calculation, the aircraft needs to be transformed from the geodetic coordinates to the body coordinates. According to the coordinate transformation theorem and Euler angle theorem, the body coordinate system can be determined by left-multiplying the geodetic coordinate system with the coordinate transformation matrix P

[14]

[0130]

[0131] Through coordinate transformation, the relative position vector r e is expressed as Equation (20):

[0132] Δr e = r H - r M = P H r H - P M r M (20)

[0133] P M is the transformation matrix from the body coordinate system of the target UAV to the global coordinate system; P H is the transformation matrix from the body coordinate system of the airliner to the global coordinate system. Given the relative position vector and relative velocity vector in the body coordinate system, assuming that the flight trajectories of the UAV and the large aircraft are independent of each other and both follow a one-dimensional normal distribution, then the relative position vector △re also follows a Gaussian distribution

[15] , that is, Equation (21):

[0134]

[0135] When the UAV is in relative motion with the aircraft, its collision probability can be regarded as the integral ratio of a probability density function that follows a three-dimensional Gaussian distribution within the collision box area. Let D be the motion area of the UAV within the collision box, then D satisfies formula (22):

[0136]

[0137] In the formula, p is used to represent the three-dimensional space vector within the collision sphere of the passenger aircraft

[16]

[0138] Suppose the particle UAV is located at A1 at time T1 within D. At this time, the relative position vector between the manned aircraft and the UAV is r1. The ellipsoidal surface that takes the center of gravity of the manned aircraft as the center of the sphere, passes through point A1 and is three-dimensionally similar to the ellipsoidal collision box surface is defined as the collision probability surface S1 at this moment. The numerical value of the surface integral of the probability density function over this collision probability surface is defined as PT1, that is, formula (23):

[0139]

[0140] At time T n is located at A n At this time, the relative position vector between the manned aircraft and the UAV is r n , the spherical surface with the center of gravity of the manned aircraft as the center of the circle and |rn| as the radius is defined as the collision surface S n At this time, P Tn is determined by formula (24):

[0141]

[0142] If the UAV has experienced a total of T1, T2... Tn points within T, then the collision surface P at this time is expressed as formula (25):

[0143]

[0144] Since the UAV is a continuous change process within the collision box, so T is continuous. Therefore, when the UAV's position changes from A1 to T within the collision box during the experience of T n Its collision volume is transformed into a triple integral, that is, formula (26):

[0145]

[0146] D1 is the collision surface area traversed by the UAV during the time T. Further, the collision probability can be obtained from the following formula (27):

[0147]

[0148] Among them, D M is the volume of the collision box.

[0149] The integration result of the above formula can be obtained by numerical integration. The above formula gives the calculation method of the collision probability, but does not give the inner collision surface nodes and outer collision surface nodes of the UAV moving in the collision box. Therefore, a specific value needs to be found in the obtained change interval to represent the maximum collision probability. Let the current time be, and the initial position vectors of the UAV and the manned airliner be r a and r b , respectively. The relative position vector can be expressed as formula (28):

[0150] r0 = r a - r b (28)

[0151] Let the relative velocity of the UAV and the airliner be V r . Since the velocity can be regarded as constant at the moment of approaching the collision area, the relative position vector of the two aircraft at this time can be expressed as formula (29):

[0152] r (t) = r0 + v r (t) (29)

[0153] In the formula, r0 is the relative position vector of the two aircraft. Let ρ(t) = r 2 (t). To find its extreme value, its derivative function can be set to 0, that is, formula (30):

[0154]

[0155] Let the solution of the above formula be t ume . When t = t ume , the position vector must take the minimum value. At this time, the double integral within the defined area between the position of the UAV and the outer boundary of the collision box is the maximum collision probability

[18]

[19] .

[0156]

[0157] The minimum value of the relative position vector can be calculated through formulas (31) and (32), and then the maximum collision probability between the UAV and the manned aircraft at a certain relative position can be calculated through formula (27).

[0158] Example 2

[0159] 4 Simulation verification

[0160] 4.1 Collision probability and safety interval

[0161] For this research, the simulated drone uses TWIN-TAI LED with a maximum flight speed of 300 km / h, and the manned aircraft uses Boeing 737-800 with a maximum speed of 880 km / h. The flight speeds are both taken as 70% of their maximum flight speeds, and the flight altitude is 7000 m. According to the flight data of the aircraft, MATLAB is used for simulation. First, the initial longitudinal interval is set at 7 km, and then, with a step size of 1 m, it is calculated 10,000 times left and right to obtain a scatter plot of the interval and the collision probability. Finally, the scatter plot is connected by interpolation method, and finally, a relationship diagram between the collision probability and the longitudinal interval is obtained. Referring to the article "Calibration Method of Safety Interval for Multi-rotor UAVs in Free Airspace of UAVs", this research defines a collision probability of 1.5% as the acceptable minimum safety interval.

[0163] Through the interpolation method, the minimum longitudinal safety interval can be obtained as 7.2983 km. In the figure, in order to determine the reasonable interval of the safety interval, a maximum safety interval needs to be found. Assuming this maximum safety interval is Max, then the safety interval range between the drone and the manned aircraft is set as [7.2983 km, Max]. Assuming that at interval M, the collision probability between the drone and the large aircraft is P(M). To obtain the specific value of Max, a interval value can be set first, and then P(M + d), P(M + 2d), P(M + 3d), …, P(M + nd) can be iteratively calculated, where the value of d is determined by the drone speed, and n is determined by d and the relative speed of the aircraft. If when d undergoes n iterations, the value of P(M + nd) changes significantly compared with other P values, then the iteration is terminated and M + nd at this time is regarded as the Max value. Usually, the rate of change is generally characterized by the slope, but the collision probability - interval in this article is obtained by the interpolation method. Using the slope to characterize the rate of change will produce errors and the results are not easy to represent. Therefore, this article considers using the area change between [P(M + i), P(M + i + 1)] and the interval axis (horizontal axis) to characterize the rate of change. According to the flight parameters of TWIN-TAI LED and Boeing 737-800 in this research, the longitudinal interval d1 = 500 m and the vertical interval d2 = 100 m are taken, and the obtained results are as Figure 6 .

[0166] When the area is depicted as a scatter plot, it can be found that the oscillation is obvious when the longitudinal interval is in the range of [7.6982 km, 8.1937 km]. In order to further determine the oscillation stable point, let d = d / 2 for iteration. After 4 consecutive iterations, it can be found that when the longitudinal interval is in the range of [7.6982 km, 7.7292 km], the oscillation begins to become gentle. When the longitudinal interval is less than 7.7292 km, the oscillation becomes gentle, and its acceptable error is within the range of 31.25 m.

[0167] It can be considered that the maximum value of the longitudinal interval within an error range of 31.25 m is 7.7292 km.

[0170] In order to obtain the vertical interval, by changing the parameters and setting the initial iteration distance to 100 m and continuously iterating, it can be found that the collision probability level oscillation begins to become gentle at the fourth time. Within a control error range of 12.25 m, the reference value of the vertical interval is considered to be [1011.6 m, 1323.12 m].

[0173] Meanwhile, the content not described in detail in this specification belongs to the prior art well-known to those skilled in the art.

[0174] It should be noted that in this text, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements not only includes those elements but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or device.

[0175] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

[0176] The embodiments of the present application have been described above. The above description is exemplary and not exhaustive, and is not limited to the disclosed embodiments. Many modifications and changes are obvious to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The choice of terms used herein is intended to best explain the principles of the embodiments, practical applications or improvements to the technology in the market, or to enable other ordinary skill in the art to understand the disclosed embodiments.

Claims

1. A method for calculating the safety interval between a UAV and a manned aircraft, characterized in that: The following steps are involved: Step 1: First, the length, wingspan, and height of the first target are obtained, and then the length, wingspan, and height of the second target are obtained. At this time, the first target is marked as a mass point, and the error caused by the volume loss of the first target as a mass point is compensated. A collision model is established according to the first target and the second target, wherein the first target is specifically a drone and the second target is specifically a manned aircraft; Step 2: When GPS is used for positioning and navigation, the multipath effect and IMU drift generated during the communication between the first target and the second target are obtained, and then the trajectory of the velocity component changing with time is obtained. When the velocity in each direction is obtained, the actual velocity vector can be obtained by vector synthesis method; Step 3: Add a three-dimensional Gaussian distribution error in the direction of the first target velocity vector, and add a Laplace distribution error in the direction of the second target velocity vector. Then, the first target is regarded as a particle approaching the collision box in three-dimensional space. The volume of the first target in the collision box is calculated by multiple integral to represent the collision probability, and the relationship between the interval and the collision probability is obtained. Step 4: Iterate the obtained safety interval and collision probability graph by intervals, and find the period with the smallest oscillation interval as the safety interval.

2. The method for calculating the safety interval between a UAV and a manned aircraft according to claim 1, characterized in that: The specific method of compensating the error caused by the volume loss of treating the first target as a particle in step 1 is: Create a collision box based on the size of the drone body. S i <<λ i i∈x,z j∈x,y (1) Based on the classic spherical and rounded cuboid collision boxes, this article defines the collision box parameters with the drone volume as follows: Then the ellipsoid collision box in the three-dimensional coordinate system is determined by formula (5): With the center of gravity of the passenger plane as the center of mass and the ellipsoid with parameter abh as the collision box, when flying in the fused airspace, when the UAV assigned to the particle flies at a certain velocity vector and the large aircraft assigned to the collision box flies at a certain velocity vector and approaches each other, it is considered that the risk of collision occurs. When flying in the joint airspace at T, the motion trajectory of the UAV is ΩT, and the collision box area attached to the large aircraft is ΦT. The collision probability is determined by the intersection of ΩT and ΦT: C=Φ T ∩Oh T (6)。 3. The method for calculating the safety interval between a UAV and a manned aircraft according to claim 1, characterized in that: In step 2, the trajectory of the velocity component changing with time is obtained. When the velocity in each direction is obtained, the actual velocity vector can be obtained by vector synthesis. The specific method is: In the positioning principle, it can be known that when the GPS signal is superimposed with multiple effects such as shadows, its positioning error is transformed into a Gaussian process. The error ε generated by the drone in the direction of running speed follows a Gaussian normal distribution with a mean u,x of 0 and a standard deviation δ = (4*16.12)*0.5 = 0.322 km [9], that is, formula (7): For large aircraft, the main source of positioning error is navigation accuracy and meteorological conditions. According to the required performance navigation principle, there is a 95% probability that a large aircraft will deviate from the planned route within a range of n nautical miles during flight. The deviation f(∈) obeys the Laplace distribution with an expected value of 0 [10], that is, formula (8): In order to eliminate the influence of positioning error, this study adds Gaussian distribution to the UAV velocity vector and Laplace distribution to the direction of the manned aircraft velocity vector for simulation calculation. When the velocity vector of an object is v and the position along the velocity direction obeys Laplace, a stochastic differential equation can be used to describe the change of the velocity component. Assume that the velocity vector of the object is: v = (vx, vy, vz) where vx, vy, vz represent the velocity components on the xyz axis respectively; consider the differential equation of formula (9): dv=-αvdt+δdW (9) dv represents a small change in the velocity component, α is the attenuation coefficient, v is the velocity component, dt is the time element, δ is the diffusion coefficient, and dW is the differential of the Wiener process (or Brownian motion). This equation contains two terms: the attenuation term and the diffusion term. The attenuation term (-ɑv) makes the velocity component tend to a stable mean. The diffusion term (δdW) introduces random noise, causing the velocity component to fluctuate around its mean. dW is the differential of the Wiener process, representing random changes. The solution of this stochastic differential equation can be simulated by the Euler numerical method, using a discretized time step Δt to approximate the differential equation to simulate the change of the velocity component over time and obtain the trajectory of the velocity component. Assuming that the velocity vector v is expressed as v = (vx, vy, vz), these velocity components obey Gaussian distribution along the velocity direction, then the change of velocity components can be expressed as formulas (10)-(12): dv x =-α x v x dt+δ x dW x (10) dv y =-α y v y dt+δ y dW y (11) dv z =-α z v z dt+δ z dW z (12) dvx, dvy and dvz represent small changes in the velocity components on the x, y and z axes respectively, ɑx, ɑy and ɑz are the corresponding attenuation coefficients, vx, vy and vz are the velocity components, δx, δy and δz are the corresponding diffusion coefficients, dWx, dWy and dWz are the differentials of the Wiener process. This formula includes the changes in the velocity components in three dimensions. Each velocity component is affected by the attenuation term and the diffusion term. The attenuation term makes the velocity component tend to a stable mean, and the diffusion term introduces random noise, causing the velocity component to fluctuate around its mean.

4. The method for calculating the safety interval between a UAV and a manned aircraft according to claim 1, characterized in that: In the step 3, the volume of the first target in the collision box is calculated by multiple integral to represent the collision probability. The specific method of obtaining the relationship between the interval and the collision probability is: Assume that the position vectors of the UAV and manned aircraft in the fusion airspace are rHrM and the velocity vectors are Vm and Vh respectively. Because the flight speed of the UAV is relatively small compared to the flight speed of the large aircraft, in order to facilitate subsequent calculations, the process of the UAV approaching the large aircraft is regarded as a process of flying at a constant speed v. Due to the speed error of the UAV, this speed cannot be simply set as the theoretical average speed of the UAV. If the steady flight speed when the aircraft approaches is v, according to the speed obstacle avoidance principle, when the speed of the UAV changes from v to v1 within tf, formula (13) is satisfied: amax is the maximum acceleration of the composite UAV’s own characteristics. Affected by technical limitations and the frequency of various crosswinds in the fusion airspace, amax satisfies a one-dimensional Gaussian distribution with a mean of amax and a variance of 0.01: If the maximum speed of the UAV during the overall stable flight stage is Vmax and the minimum speed is Vmin, then the stable flight speed Vau of the UAV considering the error can be obtained by formula (15): The relative position vector Δr and relative velocity vector Δv between the UAV and the passenger plane at time T can be expressed as equations (16) and (17) respectively: Δr e =r H -r M (16) Δv e =v H -v M (17) The relative nominal position of the UAV and the passenger plane is given by equation (18). To perform probability calculation, the aircraft needs to be transformed from the geodetic coordinates to the body coordinates. According to the coordinate system transformation theorem and the Euler angle theorem, the body coordinate system can be determined by multiplying the coordinate transformation matrix P by the geodetic coordinate system: Through coordinate system transformation, the relative position vector re is expressed as formula (20): Δr e =r H -r M =P H r H -P M r M (20) PM is the transformation matrix from the target UAV body coordinate system to the global coordinate system; PH is the transformation matrix from the passenger aircraft body coordinate system to the global coordinate system. When the relative position vector and relative velocity vector of the body coordinate system are known, assuming that the flight trajectories of the UAV and the large aircraft are independent of each other and both obey a one-dimensional normal distribution, then the relative position vector △re also obeys a Gaussian distribution, that is, formula (21): Then, when the UAV and the aircraft move relative to each other, their collision probability can be regarded as the integral ratio of a probability density function that obeys a three-dimensional Gaussian distribution in the collision box area. Let D be the movement area of ​​the UAV in the collision box, then D satisfies formula (22): In the formula, p is used to represent the three-dimensional space vector within the collision sphere of the passenger plane; Assume that the particle drone is located at A1 at time T1 in D. At this time, the relative position vector of the manned aircraft and the drone is r1. The ellipsoid surface with the center of gravity of the manned aircraft as the sphere center passing through point A1 and similar to the ellipsoid collision box surface in three dimensions is defined as the collision probability surface S1 at this moment. The integral value of the probability density function on the surface of this collision probability surface is roughly defined as PT1, that is, formula (23): At time Tn, the aircraft is located at An. At this time, the relative position vector of the manned aircraft and the unmanned aircraft is rn. The sphere with the center of gravity of the manned aircraft as the center and |rn| as the radius is defined as the collision surface Sn at this moment. At this time, PTn is determined by formula (24): If the UAV experiences a total of T1, T2, ..., Tn points within T time, then the collision surface P at this time is expressed as formula (25): Since the drone is in a continuous changing process in the collision box, T is continuous. Therefore, when the drone changes its position from A1 to Tn in the collision box after T, its collision volume is converted into a triple integral, that is, formula (26): D1 is the collision surface area traversed by the UAV in duration T. The further collision probability can be obtained by the following formula (27): Where DM is the collision box volume; The integral result of the above formula can be obtained by numerical integration. The above formula gives the calculation method of the collision probability, but does not give the inner collision surface nodes and outer collision surface nodes of the UAV moving in the collision box. Therefore, a specific value needs to be found within the obtained variation range to represent the maximum collision probability. Assuming that the current time is, the initial position vectors of the UAV and the manned aircraft are ra and rb respectively, and the relative position vector can be expressed as formula (28): r0＝r a -rb (28) Assume that the relative speed between the UAV and the passenger plane is Vr. Since the speed can be considered constant when approaching the collision area, the relative position vectors of the two aircraft can be expressed as formula (29): r (t) =r0+v r (t) (29) Where r0 is the relative position vector of the two machines. Let ρ(t) = r2(t). To obtain its extreme value, its derivative function can be set to 0, that is, formula (30): Assume that the solution of the above formula is tume. When t = tume, the position vector must reach the minimum value. At this time, the multi-integral of the position of the drone and the area defined by the outer boundary of the collision box is the maximum collision probability. The minimum value of the relative position vector can be calculated by formula (31) (32), and then the maximum collision probability between the UAV and the manned aircraft at a certain relative position can be calculated by formula (27).

Images