An airport risk area determination method, electronic equipment and storage medium

Abstract

The application provides an airport risk area determination method, electronic equipment and storage medium, relates to the field of airport risk area determination, and the method comprises the steps of: obtaining each target historical flight trajectory point information of all aircrafts taking off and landing in a target airport in a historical time period; obtaining a sub-region set Q; obtaining a trajectory point density set lambda; obtaining a vertical distance list set D; obtaining a target vertical distance corresponding to each sub-region to obtain a target vertical distance set MD; obtaining a loss of control index SK of a target unmanned aerial vehicle; determining a collision risk index of each sub-region according to the lambda, the MD and the SK to obtain a collision risk index set R; determining each sub-region in the Q as a risk area corresponding to a risk level; setting a corresponding color for each sub-region; the risk area divided by the application can match preset flight data of the target unmanned aerial vehicle, and the user can very intuitively see the risk area with a higher risk level.

Description

Technical Field

[0001] This invention relates to the field of airport risk area determination, and in particular to a method, electronic device, and storage medium for determining airport risk areas. Background Technology

[0002] To prevent the multiple risks posed by drones operating around airports, especially collisions between drones and aircraft, geofencing has become widely used to define protective boundaries around geographic points where unmanned aerial vehicles should not intrude. In existing technologies, drone activity around airports can lead to serious safety hazards. However, existing methods for delineating airport risk zones are often subjective assumptions without sufficient data support and do not incorporate drone flight information. As a result, users flying drones around airports cannot intuitively see airport risk areas that accurately match their current drone flight information. Summary of the Invention

[0003] To address the aforementioned technical problems, the technical solution adopted by this invention is as follows:

[0004] According to a first aspect of this application, a method for determining airport risk areas is provided, the method comprising the following steps:

[0005] S100: Obtain the historical flight trajectory information of each target for all aircraft that have taken off and landed at the target airport within the historical time period; the historical flight trajectory point corresponding to any target historical flight trajectory point is located within the preset area W corresponding to the target airport.

[0006] S200, divide W into several adjacent sub-regions to obtain the sub-region set Q = (Q1, Q2, ..., Q...). a Q b ), a = 1, 2, ..., b; where Q a Let b be the a-th sub-region that divides W into several sub-regions, and b be the number of sub-regions obtained after dividing W.

[0007] S300, based on Q, obtain the trajectory point density corresponding to each sub-region to obtain the trajectory point density set λ = (λ1, λ2, ..., λ300). a , …, λ b ); where λ a For Q a The corresponding trajectory point density; λ a =S a / M a S a For Q a The number of historical flight trajectory points of internal targets, M a For Q a The area;

[0008] S400, obtain the vertical distance between each target's historical flight trajectory point and the target UAV's preset flight path RT, to obtain a vertical distance list set D = (D1, D2, ..., D...). a D b ); where D a For Q a The corresponding vertical distance list; D a =(D a，1 D a，2 D a，i D a，f(a) ), i=1, 2,..., f(a); D a，i For Q a The vertical distance between the i-th target's historical flight trajectory point and RT is f(a), where Q is the vertical distance between them. a The number of historical flight trajectory points of the internal target;

[0009] S500, based on D, obtain the target vertical distance corresponding to each sub-region to obtain the target vertical distance set MD = (MD1, MD2, ..., MD4). a , ...,MD b ); where MD a For Q a Corresponding vertical distance to the target; MD a =MIN(D a MIN() is the default function for finding the minimum value;

[0010] S600, obtain the runaway index SK = α × t of the target UAV; where α is the preset stability parameter of the target UAV; t is the planned flight time of the target UAV; 0 < α < 1;

[0011] S700, based on λ, MD, and SK, determine the collision risk index for each sub-region to obtain the collision risk index set R = (R1, R2, ..., R...). a , ..., R b ); where R a For Q a The corresponding collision risk index; R a =ω1×λ a +ω2×1 / MD a +ω3×SK;R a Used to characterize the target UAV at Q a The magnitude of the risk of a collision occurring within the confines;

[0012] S800, based on R, each sub-region in Q is identified as a risk area corresponding to the risk level;

[0013] S900 assigns a corresponding color to each sub-region based on its risk level.

[0014] According to another aspect of this application, a non-transitory computer-readable storage medium is also provided, wherein at least one instruction or at least one program is stored in the storage medium, and the at least one instruction or at least one program is loaded and executed by a processor to implement the above-described airport risk area determination method.

[0015] According to another aspect of this application, an electronic device is also provided, including a processor and the aforementioned non-transitory computer-readable storage medium.

[0016] The present invention has at least the following beneficial effects:

[0017] The airport risk area determination method of the present invention divides a preset area in the target airport into several sub-areas and obtains the trajectory point density corresponding to each sub-area; based on the preset flight data of the target UAV and the trajectory point density corresponding to each sub-area, the collision risk index of each sub-area is determined, and then the average collision risk index and the volatility of the collision risk index are determined; based on the relationship between the collision risk index corresponding to each sub-area and the average collision risk index and the volatility of the collision risk index, the preset area is divided into several risk areas, and a corresponding color is assigned to each risk area; since the division of risk areas is based on the target airport's historical flight trajectory point information and the preset flight data of the target UAV, the division of risk areas in the present invention is data-supported and can be matched with the preset flight data of the target UAV. Assigning a corresponding color to each risk area allows users to intuitively see risk areas with higher risk levels. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 A flowchart of an airport risk area determination method provided in an embodiment of the present invention;

[0020] Figure 2 The decision framework diagram corresponding to the AHP (Analytic Hierarchy Process) method provided in the embodiments of the present invention is shown. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] It should be noted that various aspects of embodiments within the scope of the appended claims are described below. It will be apparent that the aspects described herein can be embodied in a wide variety of forms, and any particular structure and / or function described herein is merely illustrative. Based on this disclosure, those skilled in the art will understand that one aspect described herein can be implemented independently of any other aspect, and two or more of these aspects can be combined in various ways. For example, any number of aspects set forth herein can be used to implement the device and / or practice the method. Additionally, this device and / or method can be implemented using structures and / or functionalities other than one or more of the aspects set forth herein.

[0023] The following will refer to Figure 1 The flowchart of the airport risk area determination method introduces a method for determining airport risk areas.

[0024] The method includes the following steps:

[0025] S100: Obtain the historical flight trajectory information of each target for all aircraft that have taken off and landed at the target airport within the historical time period; the historical flight trajectory point corresponding to any target historical flight trajectory point is located within the preset area W corresponding to the target airport.

[0026] In this embodiment, for an airport, the historical flight trajectory points of aircraft are most concentrated within a certain range of the runway and its surrounding area. Therefore, the preset area W is the area within a certain range of the runway and its surrounding area. For example, the preset area W is the area outside the airport obstacle restriction surface, 10 kilometers on each side of the airport runway centerline, and 20 kilometers outside the runway end.

[0027] In this embodiment, the target's historical flight trajectory point information can be obtained through the following steps:

[0028] S110, Obtain the initial historical flight trajectory information of all aircraft taking off and landing at the target airport within the historical time period, so as to obtain the initial historical flight trajectory information set A = (A1, A2, ..., A...) corresponding to the target airport. c A d ), c = 1, 2, ..., d; where A cA represents the c-th initial historical flight trajectory information corresponding to the target area within the historical time period, and d represents the number of initial historical flight trajectory information corresponding to all aircraft taking off and landing at the target airport within the historical time period; c =(A c，1 A c，2 A c，j A c，g(c) ), j=1, 2,..., g(c); A c，j For A c Information on the j-th initial historical flight trajectory point within the range, g(c) is A c The number of initial historical flight trajectory points within the region.

[0029] In this embodiment, the aircraft can send real-time flight information to the ground control center through the ADS-B system, and each piece of flight information is recorded. Each piece of flight information is the initial historical flight trajectory point information. By obtaining the historical flight trajectory point information corresponding to each flight trajectory information of all aircraft taking off and landing at the target airport within the historical time period recorded by the ground control center, the initial historical flight trajectory information set A corresponding to the target airport can be obtained. The historical time period can be several years in the past. ADS-B is an advanced avionics communication technology. Through the ADS-B transmitter installed on the aircraft, the aircraft's status data is broadcast every second. This data includes position coordinates, altitude (altitude), speed, heading, and identification information.

[0030] S120, perform noise filtering on A to obtain the candidate historical flight trajectory information set TA = (TA1, TA2, ..., TA3) corresponding to the target airport. c , ...,TA d ); where TA c For A c Candidate historical flight trajectory information obtained after noise filtering; TA c =(TA) c，1 TA c，2 , ...,TA c，p , ...,TA c，x(c) ), p=1, 2,...,y(c); TA c，p For TA c Information on the p-th candidate historical flight trajectory point within the range, where x(c) is the TA. c The number of candidate historical flight trajectory points.

[0031] In this embodiment, there will be some noise in A, such as sensor noise, signal interference or noise generated by other unstable factors. It should be noted that those skilled in the art can use existing filtering algorithms as needed, such as using an average filtering algorithm to filter out all the noise in A to obtain the candidate historical flight trajectory information set TA corresponding to the target airport, which will not be elaborated here.

[0032] S130, iterate through TA, if TA c，p The distance D between the corresponding candidate historical flight trajectory point and the preset center point within the target airport c，p <D0, then TA c，p The intermediate historical flight trajectory points are identified to obtain the intermediate historical flight trajectory information set ZA = (ZA1, ZA2, ..., ZA...). c , ..., ZA d ); among which, ZA c For TA c Corresponding intermediate historical flight trajectory information; ZA c =(ZA c，1 ZA c，2 , ..., ZA c，q , ..., ZA c，y(c) ), q=1, 2,…,y(c); ZA c，q for ZA c Information on the q-th intermediate historical flight trajectory point within the range, y(c) is ZA c The number of historical flight trajectory points within the inner middle; where D0 is a preset distance threshold; D c，p According to TA c，p The longitude and latitude of the point, as well as the longitude and latitude of the preset center point, are obtained.

[0033] In this embodiment, the preset center point within the target airport can be the airport's control tower, D. c，p It can be obtained through the following steps:

[0034] S131, TA c，p longitude lon_c p and latitude lat_c p And convert the longitude lon0 and latitude lat0 of the preset center point within the target airport into radians.

[0035] lon_radc p =lon_c p ×π / 180; lat_radc p =lat_c p ×π / 180; lon_rad0=lon0×π / 180; lat_rad0=lat0

[0036] ×π / 180; where lon_radc p for lon_c p The radian representation, lat_radc p for lat_c p The radian representation is given by lon_rad0, where lon_rad0 is the radian representation of lon0, and lat_rad0 is the radian representation of lat0.

[0037] S132, Calculate TA using the Haversine formula. c，p The distance between the target airport and the preset center point.

[0038] TA c，p The distance D between the target airport and the preset center point c，p = 2R×arcsin[sin 2 ((lat_rad0-lat_radc p ) / 2)+cos(lat_rad0)×cos(lat_radc p )×sin 2 ((lon_rad0-lon_radc p ) / 2)] 1 / 2 Where R is the Earth's average radius (approximately 6371 km).

[0039] S140, Traverse ZA, and use a preset interpolation algorithm to perform interpolation processing on the adjacent intermediate historical flight trajectory points in each intermediate historical flight trajectory information to obtain the interpolation point information list set EA = (EA1, EA2, ..., EA) corresponding to ZA. c , ..., EA d ); where EA c for ZA c The corresponding list of interpolation points; EA c =(EA) c，1 EA c，2 , ..., EA c，r , ..., EA c，y(c)-1 ), r=1, 2,..., y(c)-1; EA c，r for ZA c，r with ZA c，r+1 Interpolation point information between; EA c，r =(Ec r，lat Ec r，lon Ec r，h ), Ec r，lat For EA c，r The latitude of the corresponding interpolation point, Ec r，lon For EA c，r The longitude of the corresponding interpolation point, Ec r，h For EAc，r The height of the corresponding interpolation point; Ec r，lat =Zc r，lat +1 / (tc r+1 -tc r )×(Zc r+1，lat -Zc r，lat )×(Mtc r -tc r ), Ec r，lon =Zc r，lon +1 / (tc r+1 -tc r )×(Zc r+1，lon -Zc r，lon )×(Mtc r -tc r ), Ec r，h =Ec r，h =Zc r，h +1 / (tc r+1 -tc r )×(Zc r+1，h -Zc r，h )×(Mtc r -tc r );Zc r，lat for ZA c，r The latitude of the corresponding historical flight trajectory point, Zc r，lon for ZA c，r The longitude of the corresponding historical flight trajectory point of the target, Zc r，h ZAc r The altitude of the corresponding historical flight trajectory point of the target, tc r for ZA c，r The time point of generation, Mtc r For tc r and tc r+1 The midpoint of time.

[0040] S150, each interpolation point information in EA and each intermediate historical flight trajectory point information in ZA are determined as the target historical flight trajectory point information.

[0041] In this embodiment, the aircraft's speed changes significantly during takeoff and landing, and ADS-B broadcasts flight trajectory information once per second, causing ZA... c The continuity of historical flight trajectory point information for the target is poor; therefore, in this embodiment, a linear interpolation algorithm is used to add ZA. c，r with ZA c，r+1 Interpolation points between them are used to improve the continuity of historical flight trajectory information of the target.

[0042] S200, divide W into several adjacent sub-regions to obtain the sub-region set Q = (Q1, Q2, ..., Q...). a Q b ), a = 1, 2, ..., b; where Q a Let b be the a-th sub-region that divides W into several sub-regions, and b be the number of sub-regions obtained after dividing W.

[0043] In this embodiment, W can be divided into grids, with each grid corresponding to a sub-region; specifically, horizontal and vertical lines with equal spacing can be added to W to obtain several sub-regions.

[0044] S300, based on Q, obtain the trajectory point density corresponding to each sub-region to obtain the trajectory point density set λ = (λ1, λ2, ..., λ300). a , …, λ b ); where λ a For Q a The corresponding trajectory point density; λ a =S a / M a S a For Q a The number of historical flight trajectory points of internal targets, M a For Q a The area.

[0045] For each sub-region, the actual area of ​​each sub-region can be obtained based on the map scale and the size of each sub-region on the map; the number of historical flight trajectory points of the target within each sub-region can be obtained by comparing the coordinates of each historical flight trajectory point of the target with the boundary coordinates of each sub-region.

[0046] S400, obtain the vertical distance between each target's historical flight trajectory point and the target UAV's preset flight path RT, to obtain a vertical distance list set D = (D1, D2, ..., D...). a D b ); where D a For Q a The corresponding vertical distance list; D a =(D a，1 D a，2 D a，i D a，f(a) ), i=1, 2,..., f(a); D a，i For Q a The vertical distance between the i-th target's historical flight trajectory point and RT is f(a), where Q is the vertical distance between them. a The number of historical flight trajectory points of the internal target.

[0047] In this embodiment, the flight path information of the target drone can be obtained, thereby determining the flight path of the target drone. Each sub-region includes several historical flight trajectory points of the target, and Q is obtained. a The vertical distance between each target's historical trajectory point and the preset flight path RT can be used to obtain D. a .

[0048] S500, based on D, obtain the target vertical distance corresponding to each sub-region to obtain the target vertical distance set MD = (MD1, MD2, ..., MD4). a , ...,MD b ); where MD a For Q a Corresponding vertical distance to the target; MD a =MIN(D a ); MIN() is the default function for finding the minimum value.

[0049] In this embodiment, Q a The minimum vertical distance among the vertical distances between each target's historical trajectory point and the preset flight path RT is determined as MD. a .

[0050] S600, obtain the runaway index SK = α × t of the target UAV; where α is the preset stability parameter of the target UAV; t is the planned flight time of the target UAV; 0 < α < 1.

[0051] In this embodiment, the performance data of the target drone can be obtained based on the flight plan data of the target drone. The performance data of the target drone includes basic information such as the weight of the target drone, the system stability parameter λ of the target drone, the maximum flight altitude of the target drone, the maximum operating speed of the target drone, and the battery life of the target drone.

[0052] It is understandable that the longer the planned flight time t, the greater the runaway index p of the target drone.

[0053] S700, based on λ, MD, and SK, determine the collision risk index for each sub-region to obtain the collision risk index set R = (R1, R2, ..., R...). a , ..., R b ); where R a For Q a The corresponding collision risk index; R a =ω1×λ a +ω2×1 / MD a +ω3×SK;R a Used to characterize the target UAV at Q a The risk of a collision occurring within the space.

[0054] In this embodiment, when determining Q a The corresponding collision risk index R a At that time, Q was taken into account. a The minimum vertical distance Q among all vertical distances between each target's historical flight trajectory point and the preset flight path information corresponding to the preset flight path is... a The corresponding trajectory point density and the runaway index of the target UAV make the determined R a This is more reasonable, and the collision risk index of each sub-region can be matched with the planned flight information of the target drone.

[0055] In this embodiment, ω1, ω2, and ω3 are determined through the following steps:

[0056] S710 uses the AHP (Analytic Hierarchy Process) to establish a hierarchical structure diagram; the hierarchical structure diagram includes, from top to bottom, the target layer, the criterion layer, and the scheme layer; the target layer is the weight ratio of the collision risk index, and the criterion layer is the trajectory point density, the reciprocal of the vertical distance to the target, and the runaway index.

[0057] AHP (Analytical Hierarchy Process): This method combines quantitative and qualitative analysis, uses the decision-maker's experience to judge the relative importance of the standards between the achievable goals, and reasonably gives the weight of each standard for each decision option, using the weights to determine the superiority or inferiority of each option.

[0058] The S720 uses a 1-9 scale method to determine the scores for trajectory point density, the reciprocal of the target vertical distance, and the runaway index.

[0059] See Figure 2 The target layer is the weight ratio of the collision risk index; the criterion layer is the trajectory point density, the reciprocal of the vertical distance to the target, and the runaway index; the scheme layer is the weight of the trajectory point density, the weight of the reciprocal of the vertical distance to the target, and the weight of the runaway index; a 1-9 scale method is established as shown in Table 1.

[0060] Table 1 1-9 Scale Method

[0061]

[0062] By combining expert scoring with three selection criteria, we aim to determine the corresponding weights of the three factors to be considered; the expert scoring results shown in Table 2 are obtained.

[0063] Table 2 Expert Scoring Results

[0064]

[0065] The results show that the reciprocal of the target vertical distance is more important than the trajectory point density, so it scores 3 points; conversely, the trajectory point density scores 0.333 points compared to the reciprocal of the target vertical distance; the trajectory point density is slightly more important than the runaway index, scoring 4 points; and the reciprocal of the target vertical distance is significantly more important than the runaway index, scoring 5 points; thus, the AHP hierarchical analysis results shown in Table 3 are obtained.

[0066] S730, based on the scoring results, obtains the AHP hierarchical analysis results; among which, the AHP hierarchical analysis results include the weight ω1 corresponding to the trajectory point density, the weight ω2 corresponding to the reciprocal of the target vertical distance, the weight ω3 corresponding to the runaway index, the maximum eigenvalue, and the CI value.

[0067] Table 3 Results of AHP (Analog-Philosophical Analysis)

[0068]

[0069] Table 3 shows that a 3rd-order judgment matrix was constructed for the trajectory point density, the reciprocal of the target vertical distance, and the runaway index (calculation method: sum-product method). The resulting eigenvectors are (0.853, 1.858, 0.289), with corresponding weights of 28.423%, 61.935%, and 9.642%. Furthermore, the largest eigenvalue, calculated from the eigenvectors, is 3.087. The CI value (0.043) is then calculated using this largest eigenvalue; where CI = (largest eigenvalue - n) / (n-1). The CI value is used for consistency testing.

[0070] S740, obtain the random consistency RI table.

[0071] In this embodiment, it should be noted that the random consistency RI table is a known table and can be obtained directly; specifically, the random consistency RI table is shown in Table 4.

[0072] Table 4 Random Consistency RI Table

[0073]

[0074] S750, obtain the RI value according to the random consistency RI table.

[0075] In this embodiment, when using the AHP (Analytic Hierarchy Process) to calculate weights, a consistency check analysis is required. The consistency check requires the use of two indicators, CI and RI. The CI value has already been calculated, and the RI value can be obtained by looking up Table 4. Thus, the consistency check results are summarized in Table 5.

[0076] S760, based on the maximum eigenvalue, CI value, and RI value, obtains the CR value.

[0077] Table 5 Summary of Consistency Test Results

[0078]

[0079] Generally, the smaller the CR value, the better the consistency of the judgment matrix. In general, if the CR value is less than 0.1, the judgment matrix satisfies the consistency test; if the CR value is greater than 0.1, it indicates that there is no consistency, and the judgment matrix should be adjusted appropriately before re-analysis.

[0080] S770, if the CR value is less than 0.1, then the weight corresponding to the trajectory point density is determined as ω1, the weight corresponding to the reciprocal of the vertical distance to the target is determined as ω2, and the weight corresponding to the runaway index is determined as ω3.

[0081] In this embodiment, the CI value calculated for the 3rd order judgment matrix is ​​0.043, and the RI value obtained by querying Table 4 is 0.520. Therefore, the calculated CR value is CI / RI = 0.083 < 0.1. Thus, the judgment matrix satisfies the consistency test, and the calculated weights are consistent. That is, the final weight values ​​corresponding to the trajectory point density are the first preset weights ω1 = 28.423%, the weight values ​​corresponding to the reciprocal of the target vertical distance are the second preset weights ω2 = 61.935%, and the weight values ​​corresponding to the runaway index are the third preset weights ω3 = 9.642%. When ω1 = 28.423%, ω2 = 61.935%, and ω3 = 9.642%, the collision risk index is relatively accurate.

[0082] S800, based on R, each sub-region in Q is identified as a risk area corresponding to the risk level.

[0083] In this embodiment, each sub-region in Q is determined as a risk region corresponding to the risk level based on R, including the following steps:

[0084] S810, Based on R, determine the mean of the average collision risk index corresponding to R and the volatility std of the collision risk index in R; where, mean = 1 / b × ∑ b a＝1 R a std = [1 / b × ∑ b a＝1 (R a -mean) 2 ] 1 / 2 .

[0085] S820, iterate through R, if mean - std ≤ R a If ≤mean, then Qa It has been identified as a Category 1 risk area.

[0086] S830, if mean < R a If Q < mean + std, then Q a It has been identified as a second-type risk area.

[0087] S840, if R a If Q ≥ mean + std, then Q a The area is identified as a third-type risk area; where η3 > η2 > η1; η1 is the risk level of the first-type risk area, η2 is the risk level of the second-type risk area, and η3 is the risk level of the third-type risk area.

[0088] S900 assigns a corresponding color to each sub-region based on its risk level.

[0089] In this embodiment, setting a corresponding color for each sub-region includes the following steps:

[0090] S910, iterate through Q, if Q a If it is a type 1 risk area, then Q will be... a Set as the first preset color;

[0091] S920, if Q a If it is a second-type risk area, then Q will be... a Set to the second preset color;

[0092] S930, if Q a If it is a third-type risk area, then Q will be... a Set to the third preset color.

[0093] In this embodiment, the first preset color can be blue, the second preset color can be yellow, and the third preset color can be red. Red risk zones prohibit drone takeoff, yellow risk zones restrict altitude, and blue zones have no restrictions. It should be noted that if R... a If <mean-std, then Q a It has been identified as a non-risk area, and Q has been designated as such. a Set it to the fourth preset color, for example, set Q a Set to light blue.

[0094] Furthermore, after step S100 and before step S200, the method further includes the following steps:

[0095] S110, obtain the density threshold MinPts; where the initial value of MinPts is from 1%*NUM to 5%*NUM, and NUM is the number of historical flight trajectory points of the target.

[0096] S120, randomly select a historical flight trajectory point U of a target that has not been visited before.

[0097] S130, obtain the number NUM of historical flight trajectory points of the target within the ε-neighborhood of U. U .

[0098] S140, if NUM U If the value is greater than or equal to MinPts, then U is marked as a core point; otherwise, U is marked as a noise point and marked as visited.

[0099] S150, obtain all historical flight trajectory points of targets within the ε-neighborhood of U to obtain the cluster corresponding to U, and mark each historical flight trajectory point of targets within the cluster corresponding to U as visited.

[0100] S160, Repeat steps S120-S150 to obtain several clusters and obtain the DBI index corresponding to the clustering results under the current MinPts.

[0101] In this embodiment, through the above steps S110-S160, several clusters corresponding to the current MinPts can be obtained. However, the density and separation of the clusters obtained by clustering under the current MinPts may not be optimal. Therefore, in this embodiment, the Davidson-Bauertin index (DBI) is used to measure the clustering results each time.

[0102] S170, adjust the value of MinPts within the range of 1%*NUM to 5%*NUM, and obtain the DBI index corresponding to the clustering result under each MinPts.

[0103] In this embodiment, for example, the first MinPts value is 1%*NUM, which can obtain the clustering result when MinPts = 1%*NUM, and also corresponds to a DBI index; the second MinPts value is 1.2%*NUM, which can obtain the clustering result when MinPts = 1.2%*NUM, and also corresponds to a DBI index; thus, multiple DBI indices can be obtained.

[0104] Specifically, for any given clustering result, calculating the DBI index corresponding to that clustering result includes the following steps:

[0105] S171, Obtain the center point of the target cluster.

[0106] In this embodiment, the target cluster is one of all clusters, and the center point of the target cluster is the average of the coordinates corresponding to all historical flight trajectory points of the target within the target cluster.

[0107] S172, obtain the average intra-cluster distance of the target cluster.

[0108] For a target cluster, calculate the distance between the historical flight trajectory point of each target within the cluster and the center point of the cluster, and calculate the average of these distances to obtain the average intra-cluster distance of the target cluster.

[0109] S173, obtain the inter-cluster distance between the target cluster and any other cluster.

[0110] In this embodiment, the inter-cluster distance between the target cluster and any other cluster is the distance between the center points of the two clusters.

[0111] S174. Determine the DBI index of the target cluster based on the average intra-cluster distance of the target cluster and the inter-cluster distance between the target cluster and any other cluster.

[0112] S175, calculate the average of the DBI indices for all clusters to obtain the DBI index corresponding to this clustering result.

[0113] S180, determine the clustering result corresponding to the smallest DBI index among all DBI indices as the target clustering result, so as to obtain the target cluster list H = (H1, H2, ..., H... j H m ), j = 1, 2, ..., m; where, H j Let m be the j-th target cluster obtained by clustering all historical flight trajectory points of the targets, and m be the number of target clusters obtained by clustering all historical flight trajectory points of the targets.

[0114] In this embodiment, the Davidson-Bourdin index (DBI) is used as a clustering quality evaluation index. DBI is a commonly used internal evaluation index used to measure the compactness and separation of clustering results. To determine whether clustering is effective, the more compact the data points within a cluster and the more separated the data points between clusters, the smaller the DBI value and the better the clustering result.

[0115] S190, set all historical trajectory points of targets within each target cluster to display in the same color, and ensure that the colors of historical trajectory points of targets within any two target clusters are different.

[0116] The airport risk area determination method in this embodiment divides a preset area in the target airport into several sub-areas and obtains the trajectory point density corresponding to each sub-area. Based on the preset flight data of the target UAV and the trajectory point density corresponding to each sub-area, the collision risk index of each sub-area is determined, and then the average collision risk index and the volatility of the collision risk index are determined. Based on the relationship between the collision risk index corresponding to each sub-area and the average collision risk index and the volatility of the collision risk index, the preset area is divided into several risk areas, and a corresponding color is assigned to each risk area. Since the division of risk areas is based on the target airport's historical flight trajectory point information and the preset flight data of the target UAV, the division of risk areas in this embodiment is data-supported and can match the preset flight data of the target UAV. Assigning a corresponding color to each risk area allows users to intuitively see risk areas with higher risk levels.

[0117] Furthermore, this embodiment also obtains the historical flight trajectory information of all aircraft that have taken off and landed at the target airport within a historical time period, clusters all the historical flight trajectory points of the targets to obtain several clusters, and adds different colors to all the historical flight trajectory points of the targets corresponding to each cluster; different colors correspond to different densities of historical flight trajectory points of the targets, so that users can first see the area with a higher density of flight trajectory points on the map, and can avoid the area when operating the drone, thereby reducing the probability of collision between the drone and the aircraft.

[0118] Furthermore, compared to the common practice of restricting the height of objects in a plane or inclined plane, this method is more targeted and provides a theoretical basis for the future integration of drones and manned aircraft in the airspace and the setting of safe intervals for electronic fences.

[0119] Furthermore, although the steps of the method in this disclosure are described in a specific order in the accompanying drawings, this does not require or imply that the steps must be performed in that specific order, or that all the steps shown must be performed to achieve the desired result. Additional or alternative steps may be omitted, multiple steps may be combined into one step, and / or a step may be broken down into multiple steps.

[0120] Embodiments of the present invention also provide a non-transitory computer-readable storage medium that can be disposed in an electronic device to store at least one instruction or at least one program related to implementing a method in the method embodiments, wherein the at least one instruction or the at least one program is loaded and executed by the processor to implement the method provided in the above embodiments.

[0121] The program product may employ any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0122] Computer-readable signal media may include data signals propagated in baseband or as part of a carrier wave, carrying readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. A readable signal medium may also be any readable medium other than a readable storage medium, capable of sending, propagating, or transmitting programs for use by or in conjunction with an instruction execution system, apparatus, or device.

[0123] The program code contained on the readable medium may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, etc., or any suitable combination thereof.

[0124] Program code for performing the operations of this application can be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, and conventional procedural programming languages ​​such as C or similar languages. The program code can execute entirely on the user's computing device, partially on the user's device, as a standalone software package, partially on the user's computing device and partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing device can be connected to the user's computing device via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computing device (e.g., via the Internet using an Internet service provider).

[0125] Embodiments of the present invention also provide an electronic device, including a processor and the aforementioned non-transitory computer-readable storage medium.

[0126] An electronic device according to this embodiment of the present application. The electronic device is merely an example and should not be construed as limiting the functionality and scope of use of the embodiments of this application.

[0127] Electronic devices are manifested in the form of general-purpose computing devices. Components of an electronic device may include, but are not limited to: at least one processor, at least one memory, and buses connecting different system components (including memory and processor).

[0128] The memory stores program code that can be executed by the processor, causing the processor to perform the steps of the various embodiments described in this specification.

[0129] The storage may include readable media in the form of volatile storage, such as random access memory (RAM) and / or cache memory, and may further include read-only memory (ROM).

[0130] The storage may also include programs / utilities having a set (at least one) of program modules, including but not limited to: an operating system, one or more applications, other program modules, and program data, each or some combination of these examples may include an implementation of a network environment.

[0131] A bus can represent one or more of several bus architectures, including a memory bus or memory controller, a peripheral bus, a graphics acceleration port, a processor, or a local bus that uses any of the various bus architectures.

[0132] The electronic device can also communicate with one or more external devices (e.g., keyboards, pointing devices, Bluetooth devices, etc.), one or more devices that enable a user to interact with the electronic device, and / or any device that enables the electronic device to communicate with one or more other computing devices (e.g., routers, modems, etc.). This communication can be performed via input / output (I / O) interfaces. Furthermore, the electronic device can communicate with one or more networks (e.g., local area networks (LANs), wide area networks (WANs), and / or public networks, such as the Internet) via a network adapter. The network adapter communicates with other modules of the electronic device via a bus. It should be understood that, although not shown in the figures, other hardware and / or software modules can be used in conjunction with the electronic device, including but not limited to: microcode, device drivers, redundant processors, external disk drive arrays, RAID systems, tape drives, and data backup storage systems.

[0133] From the above description of the embodiments, those skilled in the art will readily understand that the exemplary embodiments described herein can be implemented by software or by combining software with necessary hardware. Therefore, the technical solutions according to the embodiments of this disclosure can be embodied in the form of a software product, which can be stored in a non-volatile storage medium (such as a CD-ROM, USB flash drive, external hard drive, etc.) or on a network, including several instructions to cause a computing device (such as a personal computer, server, terminal device, or network device, etc.) to execute the methods according to the embodiments of this disclosure.

[0134] Embodiments of the present invention also provide a computer program product including program code, which, when the program product is run on an electronic device, causes the electronic device to perform the steps of the methods described above in various exemplary embodiments of the present invention.

[0135] While specific embodiments of the invention have been described in detail by way of example, those skilled in the art should understand that the above examples are for illustrative purposes only and are not intended to limit the scope of the invention. Those skilled in the art should also understand that various modifications can be made to the embodiments without departing from the scope and spirit of the invention. The scope of the invention is defined by the appended claims.

Claims

1. A method for determining airport risk areas, characterized in that, The method includes the following steps: S100: Obtain the historical flight trajectory information of each target for all aircraft that have taken off and landed at the target airport within the historical time period; the historical flight trajectory point corresponding to any target historical flight trajectory point is located within the preset area W corresponding to the target airport. S200, divide W into several adjacent sub-regions to obtain a sub-region set Q = (Q1, Q2, ..., Q...). a Q b ), a = 1, 2, ..., b; where Q a Let b be the a-th sub-region that divides W into several sub-regions, and b be the number of sub-regions obtained after dividing W. S300, based on Q, obtain the trajectory point density corresponding to each sub-region to obtain the trajectory point density set λ=(λ1, λ2, ..., λ a , …, λ b ); where λ a For Q a The corresponding trajectory point density; λ a =S a / M a S a For Q a The number of historical flight trajectory points of internal targets, M a For Q a The area; S400, obtain the vertical distance between each target's historical flight trajectory point and the target UAV's preset flight path RT, to obtain a vertical distance list set D=(D1, D2, ..., D... a D b ); where D a For Q a The corresponding vertical distance list; D a =(D a，1 D a，2 D a，i D a，f(a) ), i=1,2,…,f(a); D a，i For Q a The vertical distance between the i-th target's historical flight trajectory point and RT is f(a), where Q is the vertical distance between them. a The number of historical flight trajectory points of the internal target; S500, based on D, obtain the target vertical distance corresponding to each sub-region to obtain the target vertical distance set MD=(MD1, MD2, ..., MD... a , ...,MD b ); where MD a For Q a Corresponding vertical distance to the target; MD a =MIN(D a MIN() is the default function for finding the minimum value; S600, obtain the runaway index of the target UAV SK=α×t; where α is the preset stability parameter of the target UAV; t is the planned flight time of the target UAV; 0<α<1; S700, based on λ, MD, and SK, determine the collision risk index for each sub-region to obtain the collision risk index set R = (R1, R2, ..., R...). a , ..., R b ); where R a For Q a The corresponding collision risk index; R a =ω1×λ a +ω2×1 / MD a +ω3×SK;R a Used to characterize the target UAV at Q a The risk of collision within the time frame; ω1 is the weight corresponding to the trajectory point density, ω2 is the weight corresponding to the reciprocal of the vertical distance to the target, and ω3 is the weight corresponding to the runaway index; S800, based on R, each sub-region in Q is identified as a risk area corresponding to the risk level; S900 assigns a corresponding color to each sub-region based on its risk level.

2. The airport risk area determination method according to claim 1, characterized in that, Step S800 includes the following steps: S810, Based on R, determine the mean of the average collision risk index corresponding to R and the volatility std of the collision risk index in R; where, mean = 1 / b × ∑ b a=1 R a std=[1 / b×∑ b a=1 (R a -mean) 2 ] 1 / 2 ; S820, iterate through R, if mean - std ≤ R a If ≤mean, then Q a This area has been identified as a Category 1 risk area. S830, if mean < R a If Q < mean + std, then Q a It has been identified as a type II risk area; S840, if R a If Q ≥ mean + std, then Q a The area is identified as a third-type risk area; where η3 > η2 > η1; η1 is the risk level of the first-type risk area, η2 is the risk level of the second-type risk area, and η3 is the risk level of the third-type risk area.

3. The airport risk area determination method according to claim 2, characterized in that, Step S900 includes the following steps: S910, iterate through Q, if Q a If it is a type 1 risk area, then Q will be... a Set as the first preset color; S920, if Q a If it is a second-type risk area, then Q will be... a Set to the second preset color; S930, if Q a If it is a third-type risk area, then Q will be... a Set to the third preset color.

4. The airport risk area determination method according to claim 1, characterized in that, The target's historical flight trajectory information is obtained through the following steps: S110, Obtain the initial historical flight trajectory information of all aircraft taking off and landing at the target airport within the historical time period, so as to obtain the initial historical flight trajectory information set A=(A1, A2, ..., A...) for the target airport. c A d ), c=1,2,…,d; where A c A represents the c-th initial historical flight trajectory information corresponding to the target area within the historical time period, and d represents the number of initial historical flight trajectory information corresponding to all aircraft taking off and landing at the target airport within the historical time period; c =(A c，1 A c，2 A c，j A c，g(c) ), j=1,2,…,g(c); A c，j For A c Information on the j-th initial historical flight trajectory point within the range, g(c) is A c The number of initial historical flight trajectory points within the region; S120, perform noise filtering on A to obtain the candidate historical flight trajectory information set TA=(TA1, TA2, ..., TA3) for the target airport. c , ...,TA d ); where TA c For A c Candidate historical flight trajectory information obtained after noise filtering; TA c =(TA c，1 TA c，2 , ...,TA c，p , ...,TA c，x(c) ), p=1,2,…,y(c);TA c，p For TA c Information on the p-th candidate historical flight trajectory point within the range, where x(c) is the TA. c The number of candidate historical flight trajectory points within the region; S130, iterate through TA, if TA c，p The distance D between the corresponding candidate historical flight trajectory point and the preset center point within the target airport c，p <D0, then TA c，p The intermediate historical flight trajectory points are identified to obtain the intermediate historical flight trajectory information set ZA = (ZA1, ZA2, ..., ZA...). c , ..., ZA d ); among which, ZA c For TA c Corresponding intermediate historical flight trajectory information; ZA c =(ZA c，1 ZA c，2 , ..., ZA c，q , ..., ZA c，y(c) ), q=1,2,…,y(c); ZA c，q for ZA c Information on the q-th intermediate historical flight trajectory point within the range, y(c) is ZA c The number of historical flight trajectory points within the inner middle; where D0 is a preset distance threshold; D c，p According to TA c，p The longitude and latitude of the point, as well as the longitude and latitude of the preset center point, are obtained; S140, Traverse ZA, and use a preset interpolation algorithm to perform interpolation processing on the information of adjacent intermediate historical flight trajectory points in each intermediate historical flight trajectory information to obtain the interpolation point information list set EA=(EA1, EA2, ..., EA) corresponding to ZA. c , ..., EA d ); where EA c for ZA c The corresponding list of interpolation points; EA c =(EA c，1 EA c，2 , ..., EA c，r , ..., EA c，y(c)-1 ), r=1, 2,…,y(c)-1; EA c，r for ZA c，r with ZA c，r+1 Interpolation point information between; EA c，r =(Ec r，lat Ec r，lon Ec r，h ), Ec r，lat For EA c，r The latitude of the corresponding interpolation point, Ec r，lon For EA c，r The longitude of the corresponding interpolation point, Ec r，h For EA c，r The height of the corresponding interpolation point; Ec r，lat =Zc r，lat +1 / (tc r+1 -tc r )×(Zc r+1，lat -Zc r，lat )×(Mtc r -tc r ), Ec r，lon =Zc r，lon +1 / (tc r+1 -tc r )×(Zc r+1，lon -Zc r，lon )×(Mtc r -tc r ), Ec r，h =Ec r，h =Zc r，h +1 / (tc r+1 -tc r )×(Zc r+1，h -Zc r，h )×(Mtc r -tc r );Zc r，lat for ZA c，r The latitude of the corresponding historical flight trajectory point, Zc r，lon for ZA c，r The longitude of the corresponding historical flight trajectory point of the target, Zc r，h ZAc r The altitude of the corresponding historical flight trajectory point of the target, tc r for ZA c，r The time point of generation, Mtc r For tc r and tc r+1 The midpoint of time; S150, each interpolation point information in EA and each intermediate historical flight trajectory point information in ZA are determined as the target historical flight trajectory point information.

5. The airport risk area determination method according to claim 1, characterized in that, After step S100 and before step S200, the method further includes the following steps: S110, obtain the density threshold MinPts; where the initial value of MinPts is from 1%*NUM to 5%*NUM, and NUM is the number of historical flight trajectory points of the target; S120, randomly select a historical flight trajectory point U of a target that has not been visited before; S130, obtain the number NUM of historical flight trajectory points of the target within the ε-neighborhood of U. U ; S140, if NUM U If the value is greater than or equal to MinPts, then mark U as a core point; otherwise, mark U as a noise point and mark U as visited. S150, obtain all historical flight trajectory points of targets within the ε-neighborhood of U to obtain the cluster corresponding to U, and mark each historical flight trajectory point of targets within the cluster corresponding to U as visited; S160, Repeat steps S120-S150 to obtain several clusters and obtain the DBI index corresponding to the clustering results under the current MinPts; S170, adjust the value of MinPts within the range of 1%*NUM to 5%*NUM, and obtain the DBI index corresponding to the clustering result under each MinPts; S180, determine the clustering result corresponding to the smallest DBI index among all DBI indices as the target clustering result, so as to obtain the target cluster list H=(H1, H2, ..., H j H m ), j=1,2,…,m; where, H j Let m be the j-th target cluster obtained by clustering all historical flight trajectory points of the targets, and m be the number of target clusters obtained by clustering all historical flight trajectory points of the targets. S190, set all historical trajectory points of targets within each target cluster to display in the same color, and ensure that the colors of historical trajectory points of targets within any two target clusters are different.

6. The airport risk area determination method according to claim 1, characterized in that, ω1, ω2, and ω3 are determined through the following steps: S710 uses the AHP (Analysis of Hierarchical Organization) method to establish a hierarchical structure diagram; the hierarchical structure diagram includes, from top to bottom, the target layer, the criterion layer, and the scheme layer; the target layer is the weight ratio of the collision risk index, and the criterion layer is the trajectory point density, the reciprocal of the vertical distance to the target, and the runaway index; S720 uses a 1-9 scale method to determine the scoring results of trajectory point density, the reciprocal of the target vertical distance, and the runaway index; S730, based on the scoring results, obtains the AHP hierarchical analysis results; among which, the AHP hierarchical analysis results include the weight ω1 corresponding to the trajectory point density, the weight ω2 corresponding to the reciprocal of the target vertical distance, the weight ω3 corresponding to the runaway index, the maximum eigenvalue, and the CI value; S740, obtain the random consistency RI table; S750, obtain the RI value according to the random consistency RI table; S760, the CR value is obtained based on the maximum eigenvalue, CI value, and RI value; S770, if the CR value is less than 0.1, then the weight corresponding to the trajectory point density is determined as ω1, the weight corresponding to the reciprocal of the vertical distance to the target is determined as ω2, and the weight corresponding to the runaway index is determined as ω3.

7. The airport risk area determination method according to claim 1, characterized in that, ω1=28.423%, ω2=61.935%, ω3=9.642%.

8. A non-transitory computer-readable storage medium, wherein the storage medium stores at least one instruction or at least one program segment, characterized in that, The at least one instruction or the at least one program segment is loaded and executed by the processor to implement the airport risk area determination method as described in any one of claims 1-7.

9. An electronic device, characterized in that, Includes a processor and the non-transitory computer-readable storage medium as described in claim 8.

Images