A method for clustering aircraft trajectories based on local adaptive dynamic time warping

Abstract

The present invention discloses an aircraft trajectory clustering method based on local adaptive dynamic time warping. First, the arrival and departure track data are obtained, and the quality analysis of the track data is performed; and the track data is preprocessed; secondly, the runways are distinguished, and a runway is assigned to each trajectory, and the trajectories taking off and landing on the same runway are further classified according to the corridor entrance; then a local adaptive dynamic time warping model is constructed to find the optimal path and the overall distance of the path; finally, the aircraft trajectories are clustered based on the partitioning method. The present invention draws on the idea of ​​intuitive visual distinction of trajectories, and fully amplifies the parts of the trajectory with lower similarity by artificially applying larger weights, so that the similarity of the parts with lower similarity is lower, and the influence of the similar parts on the similarity of the entire trajectory is not considered, so that the clustering results are more accurate, explanatory and more targeted.

Description

Technical Field

[0001] The present invention belongs to the technical field of automation and intelligentization of civil aviation air traffic control, and in particular relates to an aircraft trajectory clustering method based on local adaptive dynamic time warping. Background Art

[0002] In recent years, with the gradual maturity of technologies such as automatic dependent surveillance-broadcast systems and secondary radars, the acquisition of massive track data has become increasingly convenient. Analysis of a large amount of track data shows that many flights do not fly according to standard approach and departure procedures in the terminal airspace due to large flight volume, weather, military control or other reasons. This shows that standard terminal approach and departure procedures are sometimes not sufficient to define a representative route structure. Therefore, mining representative actual flight routes from a large amount of historical track data has become a hot topic in current air traffic research. At present, the most mainstream and effective way to achieve this goal is to perform cluster analysis on trajectory data. Trajectory clustering is the process of dividing a set of trajectories into several subsets with similar characteristics. Among them, the cluster center can represent the prevailing flight mode of most trajectories to a certain extent, and can provide a reference for evaluating the rationality of the formulation of standard flight procedures. By observing the representative trajectories and abnormal trajectories found by trajectory clustering, we can peek into the shortcomings of the current airspace structure, which plays a vital role in the scientific division of airspace sectors, the efficient use of airspace resources and the safety of airspace operations. In addition, track clustering is the premise and basis of track prediction, and it is also a way to generate tracks and simulate airspace operations.

[0003] How to define the similarity between objects has a significant impact on the clustering results because it determines the direction of the clustering search. Similarity is generally defined by the distance between objects. Due to the inconsistency of track lengths, there are generally two methods for processing track lengths. One is to resample, smooth and reconstruct the tracks to make them have the same length, and then use general distance formulas such as Euclidean distance to calculate the distance between tracks as similarity; the other is to not do any processing on the track length and directly use methods that can measure the similarity between tracks of different lengths, such as dynamic time warping, Fréchet distance, Hausdorff distance, etc.

[0004] However, in the clustering process, it was found that some trajectories were often classified into one category when clustering because the distances of most stages were close, but some stages were far away. Although adjusting the parameters to increase the number of categories can sometimes solve this problem, it may cause trajectories that originally belonged to one category to become two categories. Therefore, starting from the distance measurement method, the dynamic time warping method that can measure the similarity between trajectories of different lengths is improved to solve the above problem. Summary of the invention

[0005] Purpose of the invention: In view of the above-mentioned problems existing in the dynamic time warping method, the present invention proposes an aircraft trajectory clustering method based on local adaptive dynamic time warping, so that the clustering results are more accurate, explanatory and more targeted.

[0006] Technical solution: The aircraft trajectory clustering method based on local adaptive dynamic time warping described in the present invention comprises the following steps:

[0007] (1) Obtain arrival and departure track data and perform quality analysis on the track data;

[0008] (2) preprocessing the track data obtained in step (1);

[0009] (3) Differentiate the runways, assign a runway to each trajectory, and further classify the trajectories taking off and landing on the same runway according to the corridor entrances;

[0010] (4) Construct a local adaptive dynamic time warping model to find the optimal path and the overall distance of the path;

[0011] (5) Clustering of aircraft trajectories based on the partitioning method.

[0012] Furthermore, the track data described in step (1) includes overflight, take-off and landing track data.

[0013] Furthermore, the quality analysis of the track data in step (1) includes: merging all tracks of a flight within one day; adding track numbers to the flights of each day in the form of year, month, day plus track sequence number, and sorting all track points of each track in chronological order; analyzing whether there are missing, repeated and abnormal data in the data; and analyzing whether there are tracks with too few track points and too large intervals between the previous and next points.

[0014] Furthermore, the step (2) comprises the following steps:

[0015] (21) Data cleaning: For repeated track points in the data, the first repeated track point is retained and all other repeated track point data are deleted; missing values ​​in the data are deleted and interpolated;

[0016] (22) Coordinate conversion: The longitude and latitude coordinates in the secondary radar data are converted from the geographic coordinate system to the coordinate values ​​of the northeast celestial rectangular coordinate system with the airport center as the origin. east ,y north 、z up :

[0017] x east =-sinθ0×(X-x0)+cosθ0×(Y-y0) (7)

[0018]

[0019]

[0020] Where X, Y, and Z are the coordinate values ​​of the corresponding Earth-centered Earth-fixed coordinate system; (x0, y0, z0) is the Earth-centered Earth-fixed coordinate of the center of the terminal area, represents the radian corresponding to the latitude of the origin, and θ0 represents the radian corresponding to the longitude of the origin;

[0021] (23) Feature normalization: Perform maximum-minimum normalization on the cleaned track information data set and map the data into the interval [-1,1];

[0022] (25) Derived attribute construction: merge all arrival and departure tracks together; count the number of different moments; count the number of different arrival track numbers and departure track numbers at each moment, and use them as the current arrival flight volume and departure flight volume at that moment, respectively.

[0023] Furthermore, the implementation process of step (3) is as follows:

[0024] A rectangular area including the runway area and the peripheral area is drawn near each runway. A line perpendicular to the edges of the two runways is drawn at the midpoint of the runway as the dividing line of the rectangular area. The rectangular area is divided into two equal parts according to the area, and labels are added to each area to distinguish them, corresponding to each runway. The number of points in each area where all track points with heights from the runway end height to the height at the boundary of the rectangular area fall is counted respectively. The track with the most points in the area belongs to the takeoff or landing area. For the track that takes off from the same runway, the distance between its last track point and each corridor entrance is calculated. For the track that lands on the same runway, the distance between its first track point and each corridor entrance is calculated. The corridor entrance to which the track is closest belongs to.

[0025] Furthermore, the step (4) comprises the following steps:

[0026] (41) Create the initial distance matrix: Assume that trajectory T i The track point sequence set is: Track T i+1 The track point sequence set is Create an m×n path matrix where the (c th ,d th ) elements are two points and The Euclidean distance is given a weight, that is ||·||2 represents the Euclidean distance between two points; each track point All contain E coordinates x east 、N coordinate y north 、U coordinate z up , new lift rate sjl, current arrival flight volume f land 、The current number of departing flights f take , heading D and other seven attribute values;

[0027] (42) The regularized path is subject to the following constraints:

[0028] 1) Boundary point constraints: The start and end points of the regularized path must be the first and last points of the path matrix;

[0029] 2) Continuity constraint: the path can only move forward one step at a time;

[0030] 3) Monotonicity constraint: The points on the path must be monotonically distributed over time;

[0031] (43) Dynamically search for the optimal path: The regular path is L = {u1,u2,…,u G}, G represents the number of pairs of similar points, u represents the index of a pair of similar points in the regular path, represented by (c, d); the Euclidean distance sequence between similar points in the regular path is represented by D L ={w1,w2,L,w G} indicates that w i Represents the Euclidean distance between a pair of similar points;

[0032] Combined with the continuity and monotonicity constraints, the path of each grid point has only three directions: if the path has passed through the grid point (c, d), then the next grid point passed is (c+1, d), (c, d+1), or (c+1, d+1);

[0033] The path that minimizes the regularization cost:

[0034]

[0035] Define a cumulative distance, starting from point (0,0) to match two sequences, each time you reach a point, the distances calculated from all previous points will be accumulated, and after reaching the end point (m,n), this cumulative distance is the final total distance; the cumulative distance is the current grid point distance, that is, the point and The sum of the Euclidean distance of the point and the cumulative distance of the smallest neighboring element that can reach the point, thus finding the optimal path:

[0036]

[0037] (44) Local adaptive weighted similarity calculation:

[0038] After finding the optimal path, find the distances between similar points from the original distance matrix according to the path and arrange them in ascending order. Find the 40% quantile pe of the distance sequence, set the weight of the distance less than pe to 0, that is, do not consider the distance of the closer part, and set the weight of the distance greater than pe to 1.5, that is, magnify the distance of the dissimilar part in the trajectory; then, calculate the average value of the weighted distance with a weight of 1.5, which is the local adaptive DTW distance between the two trajectories.

[0039] Furthermore, the step (5) comprises the following steps:

[0040] (51) Randomly select K trajectories from the set of all trajectories as the initial representative trajectories;

[0041] (52) Assign each remaining trajectory to the nearest representative trajectory. Local adaptive dynamic time warping is used here to calculate the similarity between each remaining trajectory and each representative trajectory. Each remaining trajectory is assigned to the representative trajectory with the greatest similarity (the smallest local adaptive dynamic time warping distance);

[0042] (53) Randomly select a non-representative trajectory o random ;

[0043] (54) Calculation random Instead of representing the trajectory o j The total cost S is:

[0044]

[0045] Among them, dist(T i ,T i+1 ) represents the local adaptive dynamic time warping distance between two trajectories, are all the trajectories T and C in the dataset j The representative trajectory o j The sum of the absolute errors of are all the trajectories T and C in the dataset random The representative trajectory o random The sum of the absolute errors of

[0046] (55) If S<0, then o random Replace o j , forming a new set of K representative trajectories;

[0047] (56) Repeat steps (53)-(55) until the representative trajectory does not change;

[0048] (57) The silhouette coefficient method is used to determine the K value; the value is between -1 and 1. The closer the K value is to 1, the more reasonable the trajectory clustering is; the closer the K value is to -1, the trajectory should be clustered into other clusters; and when the K value is close to 0, the sample is on the boundary of two clusters.

[0049] Beneficial effects: Compared with the prior art, the present invention has the following beneficial effects:

[0050] 1. For trajectories with too similar trajectories or too large distances between them, traditional distance metrics may mistakenly classify them into one category, and adjusting the parameter value cannot separate them. The present invention solves the above problem by setting the distance weight of the similar parts of the trajectories to 0, that is, ignoring the similar parts, and amplifying the distance of the dissimilar parts, that is, setting a larger weight, so as to separate them.

[0051] 2. Insensitive to outliers: Since the partition-based K-center point clustering method is used, it is different from K-means clustering. The former selects the actual object as the representative trajectory, while the latter calculates the mean of all objects in the class as the representative pair. Therefore, the latter is more sensitive to outliers, while the former is not.

[0052] 3. Two features, namely the current number of incoming flights and the current number of departing flights, which have an important impact on the clustering results, are proposed. Trajectories that originally belonged to one category without this feature can be divided into multiple categories according to the current flight volume status of the airspace, making the clustering results more accurate and explanatory;

[0053] 4. Divide the trajectories into several categories according to the runways and corridor entrances, that is, separate the data set first, and then cluster each type of trajectory separately on this basis. This can make the clustering results more targeted and prevent obvious errors. For example, due to the close distance, a trajectory that does not belong to a runway (or corridor entrance) is attributed to the runway (or corridor entrance). BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 is a flow chart of the present invention;

[0055] Figure 2 This is a schematic diagram of the runway divisions;

[0056] Figure 3 A schematic diagram of the regular path. DETAILED DESCRIPTION

[0057] The present invention is further described in detail below with reference to the accompanying drawings.

[0058] The present invention proposes an aircraft trajectory clustering method based on local adaptive dynamic time warping, such as Figure 1 As shown, the specific steps include:

[0059] Step 1: Obtain arrival and departure track data and perform quality analysis on the track data.

[0060] Based on all the existing track data of a certain area, including flyover, take-off and landing tracks, the take-off and landing track data belonging to a specific terminal area are extracted, that is, all the tracks of take-off and landing at a certain airport within an area with a radius of 80 kilometers centered on the airport, including information recording time, flight number, aircraft location information (longitude, latitude and altitude), heading, ascent and descent rate, aircraft speed, etc., and the arrival and departure data are divided. All subsequent steps are for separate processing of approach and departure tracks.

[0061] Merge all the tracks belonging to a day. Considering that the take-off or landing of a flight may be distributed in the last few minutes of the previous hour and the first few minutes of the next hour, and the take-off and landing time of the flight in the terminal area varies, ranging from 15 minutes to half an hour, so before exploring the data quality, first merge the track data of 24 hours a day.

[0062] Add track numbers to distinguish all tracks. There may be multiple tracks corresponding to one flight number in a day. Therefore, for such tracks, if the time interval between two consecutive points exceeds 15 minutes, then the flight numbers are artificially distinguished, and there is one track before the interval and another track after the interval. That is, since it is impossible for an aircraft to take off immediately after landing or land immediately after taking off, for a flight number corresponding to multiple tracks, there must be a time interval greater than 15 minutes between two consecutive points. In addition, the same flight number may exist between days, months, and years. Therefore, in order to distinguish the tracks, the track numbers are added to the flights of each day in the form of year, month, and day plus the track sequence number, and all the track points of each track are sorted in chronological order. In this way, for all track data within one year, the tracks can be distinguished from each other.

[0063] Analyze whether there are missing, duplicate and abnormal data in the data. Compare all the attributes of all track points in the track data set. When all the attributes between two track points are the same, the two track points are duplicate values; check the attribute values ​​of the track points. When there are null values, the track information is missing values; clustering method is used to determine the abnormal values. Since the track information data is multivariate data, fast clustering method is used to process multivariate abnormal values, and the data objects are grouped into multiple clusters, so as to mine outliers and determine abnormal value data.

[0064] Analyze whether there are too few track points or too large intervals between the two points. Count the number of points in all tracks and the time intervals between all track points in each track, and analyze whether there are too few track points or too large intervals between the two points.

[0065] Step 2: Preprocess the track data obtained in step 1.

[0066] Data cleaning. In order to ensure the quality and efficiency of subsequent clustering, for repeated track points in the data, retain the first repeated track point and delete all other repeated track point data; the missing values ​​in the data are deleted and interpolated. When the proportion of missing attributes in the data is greater than 85%, delete the information data of the track point, otherwise use the regression method to interpolate the missing values ​​based on other non-missing attribute values; for outliers, uniformly adopt the deletion method; delete the tracks with less than 90 track points and the time interval between two consecutive points is more than 90 seconds, because the former cannot reflect the entire approach process, and the latter has too many missing track points. Supplementing by interpolation and other methods will result in excessive errors, affecting the quality of subsequent clustering;

[0067] Coordinate conversion. In order to simplify the calculation of spatial distance and present the changes of aircraft track information more intuitively, it is necessary to convert the latitude and longitude coordinates in the secondary radar data from the geographic coordinate system to the northeast celestial rectangular coordinate system with the center of the airport as the origin. First, convert the data from the geographic coordinate system to the geocentric earth-fixed rectangular coordinate system. The main coordinate change process is as follows:

[0068] 1) Calculate the ellipsoid flattening F, eccentricity e and curvature radius r corresponding to the spatial spherical coordinate system:

[0069]

[0070] e=F×(2-F) (2)

[0071]

[0072] in, represents the radian corresponding to the latitude, a represents the radius of the earth's major axis, and b represents the radius of the earth's minor axis. Then, calculate the coordinate values ​​X, Y, and Z of the corresponding Earth-centered Earth-fixed coordinate system:

[0073]

[0074]

[0075]

[0076] Among them, θ represents the radian corresponding to the longitude, and h represents the altitude.

[0077] 2) The Earth-centered Earth-fixed coordinates are converted into a Northeastern celestial rectangular coordinate system with the airport center as the origin, where the positive Z axis coincides with the normal of the ellipsoid, the positive Y axis points to the north, and the positive X axis points to the east.

[0078] The Earth-fixed coordinates (x0, y0, z0) of the center of the terminal area are taken as the origin. The radian corresponding to the latitude of the origin, θ0, the radian corresponding to the longitude of the origin, and h0, the height of the origin. Calculate the coordinate value x of the ENU rectangular coordinate system based on the coordinates of the origin east ,y north 、z up :

[0079] x east =-sinθ0×(X-x0)+cosθ0×(Y-y0) (7)

[0080]

[0081]

[0082] Feature normalization: Perform maximum-minimum normalization on the cleaned track information dataset and map the data to the [-1,1] interval.

[0083] The normalized heading, E coordinate, N coordinate, and U coordinate attributes are multiplied by 2 respectively so that each feature has an appropriate influence weight.

[0084] Derived attribute construction. Current number of approaching flights and number of departing flights: In actual operation, aircraft may be affected by the number of flights in the current airspace. When the number of flights is large, aircraft generally strictly follow the five-side approach or takeoff. When the number of flights is even larger, the controller generally implements radar guidance on the aircraft. At that time, the aircraft trajectory is more uncertain; and when the number of flights is small, the aircraft trajectory will also be uncertain. For example, in order to save unnecessary time, the controller may command the aircraft to cut directly to the last approach edge. Therefore, two attributes, the current number of approaching flights and the current number of departing flights, are constructed. The specific construction steps are as follows:

[0085] ① Combine all entry and exit trajectories together;

[0086] ② For all track data, count the number of different values ​​in the time attribute, that is, the total number of different time values ​​in the time attribute.

[0087] ③ Count the number of different arrival track numbers and departure track numbers at each moment, and use them as the current arrival and departure flight volumes at that moment.

[0088] Step 3: Differentiate the runways, assign a runway to each trajectory, and further classify the trajectories taking off and landing on the same runway according to the corridor entrances.

[0089] Runway matching. First, partition the runways. In order to determine which runway each trajectory takes off or lands on, a flat rectangular area is drawn near each runway. The rectangular area not only includes the runway area, but also part of its outer area. Draw a line perpendicular to the edges of the two runways at the midpoint of the runway as the dividing line of the rectangular area. Divide the rectangular area into two equal parts according to the area, and add labels to each area to distinguish them, corresponding to their respective runways. For example, if an airport has two runways, divide it into 4 areas and use four labels A, B, C, and D to distinguish them, such as Figure 2 As shown. Secondly, assign a runway to the trajectory. Check the approximate height of the trajectory data at the end of the runway and the boundary of the rectangular area, which are approximately 250 meters and 1000 meters respectively. Then count the number of points in each area where all track points with a height of 250 meters to 1000 meters fall. The track with the most points in the area belongs to the take-off or landing area. Track points with a height of less than 250 meters in the trajectory are not considered because the points on the runway are relatively dense when the aircraft is running at a slow speed and the radar scanning cycle remains unchanged. In addition, the points in the taxiing phase may not belong to the runway area where it takes off or lands, which will interfere with the judgment of the runway to which the trajectory belongs.

[0090] Corridor entrance matching. In order to be more targeted during clustering, the trajectories of takeoff and landing on the same runway need to be further classified according to the corridor entrances. Since only the trajectories within the 80-kilometer radius centered on the airport are selected, this area may not include some corridor entrances. Therefore, for the trajectories belonging to the same runway takeoff, calculate the distance between its last trajectory point and each corridor entrance; for the trajectories belonging to the same runway landing, calculate the distance between its first trajectory point and each corridor entrance; the corridor entrance to which the trajectory is closest belongs to.

[0091] Step 4: Build a local adaptive dynamic time warping model to find the optimal path and the overall distance of the path.

[0092] Create an initial distance matrix. Assume that trajectory T i The track point sequence set is Track T i+1 The track point sequence set is Each track point All contain E coordinates x east 、N coordinate y north 、U coordinate z up , new lift rate sjl, current arrival flight volume f land 、The current number of departing flights f take , heading D and other seven attribute values. Create an m×n path matrix. The first (c th ,d th ) elements are two points and The Euclidean distance of ||·||2 represents the Euclidean distance between two points.

[0093] Dynamically find the optimal path. The regular path is L = {u1,u2,L,u G}, G represents the number of similar points in total, and u represents the index of a pair of similar points in the regular path, represented by (c, d), such as Figure 2 The Euclidean distance sequence between similar points on the regularized path is represented by D L ={w1,w2,…,w G} indicates that w i Represents the Euclidean distance between a pair of similar points.

[0094] Regularized paths are subject to several constraints:

[0095] Boundary point constraint: The start and end points of the regularized path must be the first and last points of the path matrix.

[0096] Continuity constraint: The path can only move forward one step at a time.

[0097] Monotonicity constraint: The points on the path must be monotonically distributed over time.

[0098] Combined with the continuity and monotonicity constraints, each path to a grid point has only three directions. For example, if the path has passed through the grid point (c, d), then the next grid point passed can only be one of the following three cases: (c+1, d), (c, d+1) or (c+1, d+1).

[0099] There are exponential paths that satisfy the above constraints, but the requirement is to minimize the following regularized cost:

[0100]

[0101] Define a cumulative distance, start matching two sequences from point (0,0), and when you reach a point, the distances calculated from all previous points will be accumulated. When you reach the end point (m,n), this cumulative distance is the final total distance. The calculation method of the cumulative distance is shown in formula (12). Among them, η(c,d) represents the dynamic update distance of a grid point, Indicates the track T i Track points and track T i+1 Track points The Euclidean distance between the points is the accumulated distance between the current grid points. and The sum of the Euclidean distance of the point and the cumulative distance of the smallest neighboring element that can reach the point. Thus, the optimal path is found.

[0102]

[0103] Local adaptive weighted similarity calculation:

[0104] After finding the optimal path, find the distance between each similar point from the original distance matrix according to the path, arrange them in order from small to large, find the 40% quantile pe of the distance sequence, set the weight of the distance less than pe to 0, that is, do not consider the distance of the closer part, and set the weight of the distance greater than pe to 1.5, that is, magnify the distance of the dissimilar part in the trajectory. Then, calculate the average value of the weighted distance with a weight of 1.5, which is the local adaptive DTW distance between the two trajectories. The setting of quantiles and weights here should be based on the situation of each type of trajectory. If there are many similar parts of the trajectory, the quantile should be selected relatively large, and the setting of weights should also be based on the situation of the trajectory to be clustered. If the distance between the dissimilar parts of the trajectory to be clustered is small, the weight should be set relatively large to magnify the dissimilar parts in the trajectory, so as to obtain an ideal clustering result.

[0105] Step 5: Cluster the aircraft trajectories based on the partitioning method.

[0106] 1) Randomly select K trajectories from the set of all trajectories as the initial representative trajectories.

[0107] 2) Assign each remaining trajectory to the nearest representative trajectory. Local adaptive dynamic time warping is used here to calculate the similarity between each remaining trajectory and each representative trajectory. Each remaining trajectory is assigned to the representative trajectory with the greatest similarity (the smallest local adaptive dynamic time warping distance).

[0108] 3) Randomly select a non-representative trajectory o random .

[0109] 4) Calculation random Instead of representing the trajectory o j The total cost S is shown in formula (13).

[0110]

[0111] Among them, dist(T i ,T i+1 ) represents the local adaptive dynamic time warping distance between two trajectories, are all the trajectories T and C in the dataset j The representative trajectory o j The sum of the absolute errors of are all the trajectories T and C in the dataset random The representative trajectory o randomThe sum of the absolute errors.

[0112] 5) If S<0, then o random Replace o j , forming a new set of K representative trajectories.

[0113] 6) Repeat steps 3) to 5) until the representative trajectory does not change.

[0114] 7) Use the silhouette coefficient method to determine the K value.

[0115] The silhouette coefficient measures the quality of clustering from the perspective of intra-cluster compactness and inter-cluster separation. Its value is between -1 and 1. The closer the value is to 1, the more reasonable the trajectory clustering is; the closer the value is to -1, the more likely the trajectory should be clustered into other clusters; the closer the value is to 0, the sample is on the boundary of two clusters. The mean of all trajectory silhouette coefficients is taken as the overall silhouette coefficient, and it is used to measure the overall quality of the clustering results.

[0116] Taking the terminal area of ​​Nanjing Lukou Airport as an example, flight data from July 20, 2019 to August 11, 2019 were used, including a total of 10,311 trajectory data, including 5,073 landing trajectories and 5,238 take-off trajectories. For all trajectories of take-off and landing clusters, the silhouette coefficient method was used to compare the effects of similarity calculated using standard dynamic time warping and local adaptive dynamic time warping, which were 0.6 and 0.8 respectively. The improved dynamic time warping method - local adaptive dynamic time warping can make the overall clustering effect better than that using standard dynamic time warping.

Claims

1. A method for clustering aircraft trajectories based on local adaptive dynamic time warping, characterized in that: The following steps are involved: (1) Obtain arrival and departure track data and perform quality analysis on the track data; (2) preprocessing the track data obtained in step (1); (3) Distinguish the runways, assign a runway to each trajectory, and further classify the trajectories taking off and landing on the same runway according to the corridor entrances; (4) Construct a local adaptive dynamic time warping model to find the optimal path and the overall distance of the path; (5) Clustering of aircraft trajectories based on the partitioning method; The step (2) comprises the following steps: (21) Data cleaning: For repeated track points in the data, the first repeated track point is retained and all other repeated track point data are deleted; missing values ​​in the data are deleted and interpolated; (22) Coordinate conversion: The longitude and latitude coordinates in the secondary radar data are converted from the geographic coordinate system to the coordinate values ​​of the northeast celestial rectangular coordinate system with the airport center as the origin. east ,y north 、z up : x east =-sinθ0×(X-x0)+cosθ0×(Y-y0) (7) Where X, Y, and Z are the coordinate values ​​of the corresponding Earth-centered Earth-fixed coordinate system; (x0, y0, z0) is the Earth-centered Earth-fixed coordinate of the center of the terminal area, represents the radian corresponding to the latitude of the origin, and θ0 represents the radian corresponding to the longitude of the origin; (23) Feature normalization: Perform maximum-minimum normalization on the cleaned track information data set and map the data into the interval [-1,1]; (24) Derived attribute construction: merge all arrival and departure tracks together; count the number of different moments; count the number of different arrival track numbers and departure track numbers at each moment, and use them as the current arrival flight volume and departure flight volume at that moment; The step (4) comprises the following steps: (41) Create the initial distance matrix: Assume that the trajectory T i The track point sequence set is: Track T i+1 The set of track point sequences is Create an m×n path matrix where the (c th ,d th ) elements are two points and The Euclidean distance is given a weight, that is ||·||2 represents the Euclidean distance between two points; each track point All contain E coordinates x east 、N coordinate y north 、U coordinate z up , new lift rate sjl, current arrival flight volume f land 、The current number of departing flights f take , heading D seven attribute values; (42) The regularized path is subject to the following constraints: 1) Boundary point constraints: The start and end points of the regularized path must be the first and last points of the path matrix; 2) Continuity constraint: the path can only move forward one step at a time; 3) Monotonicity constraint: The points on the path must be monotonically distributed over time; (43) Dynamically search for the optimal path: The regular path is L = {u1,u2,…,u G }, G represents the number of pairs of similar points, u represents the index of a pair of similar points in the regular path, represented by (c, d); the Euclidean distance sequence between similar points in the regular path is represented by D L ={w1,w2,…,w G } indicates that w i Represents the Euclidean distance between a pair of similar points; Combined with the continuity and monotonicity constraints, the path of each grid point has only three directions: if the path has passed through the grid point (c, d), then the next grid point passed is (c+1, d), (c, d+1), or (c+1, d+1); The path that minimizes the regularization cost: Define a cumulative distance, starting from point (0,0) to match two sequences, each time you reach a point, the distances calculated from all previous points will be accumulated, and after reaching the end point (m,n), this cumulative distance is the final total distance; the cumulative distance is the current grid point distance, that is, the point and The sum of the Euclidean distance of the point and the cumulative distance of the smallest neighboring element that can reach the point, thus finding the optimal path: (44) Local adaptive weighted similarity calculation: After finding the optimal path, find the distances between similar points from the original distance matrix according to the path, arrange them in ascending order, find the 40% quantile pe of the distance sequence, set the weight of the distance less than pe to 0, that is, do not consider the distance of the closer part, and set the weight of the distance greater than pe to 1.5, that is, magnify the distance of the dissimilar part in the trajectory; then, calculate the average value of the weighted distance with a weight of 1.5, which is the local adaptive DTW distance between the two trajectories; The step (5) comprises the following steps: (51) Randomly select K trajectories from the set of all trajectories as the initial representative trajectories; (52) Assign each remaining trajectory to the nearest representative trajectory. Local adaptive dynamic time warping is used here to calculate the similarity between each remaining trajectory and each representative trajectory. Each remaining trajectory is assigned to the representative trajectory with the greatest similarity. (53) Randomly select a non-representative trajectory o random ; (54) Calculation random Instead of representing the trajectory o j The total cost S is: Among them, dist(T i ,T i+1 ) represents the local adaptive dynamic time warping distance between two trajectories, are all the trajectories T and C in the dataset j The representative trajectory o j The sum of the absolute errors of are all the trajectories T and C in the dataset random The representative trajectory o random The sum of the absolute errors of (55) If S < 0, then o random Replace o j , forming a new set of K representative trajectories; (56) Repeat steps (53)-(55) until the representative trajectory does not change; (57) The silhouette coefficient method is used to determine the K value; the value is between -1 and 1. The closer the K value is to 1, the more reasonable the trajectory clustering is; the closer the K value is to -1, the trajectory should be clustered into other clusters; and when the K value is close to 0, the sample is on the boundary of two clusters.

2. The aircraft trajectory clustering method based on local adaptive dynamic time warping according to claim 1, characterized in that: The track data described in step (1) includes flyover, take-off and landing track data.

3. The aircraft trajectory clustering method based on local adaptive dynamic time warping according to claim 1, characterized in that: The quality analysis of the track data described in step (1) includes: merging all tracks of a flight within one day; adding track numbers to the flights of each day in the form of year, month, day plus track sequence number, and sorting all track points of each track in chronological order; analyzing whether there are missing, repeated and abnormal data in the data; analyzing whether there are tracks with too few track points and tracks with too large intervals between the previous and next points.

4. The aircraft trajectory clustering method based on local adaptive dynamic time warping according to claim 1, characterized in that: The implementation process of step (3) is as follows: A rectangular area including the runway area and the peripheral area is drawn near each runway. A line perpendicular to the edges of the two runways is drawn at the midpoint of the runway as the dividing line of the rectangular area. The rectangular area is divided into two equal parts according to the area, and labels are added to each area to distinguish them, corresponding to each runway. The number of points in each area where all track points with heights from the runway end height to the height at the boundary of the rectangular area fall is counted respectively. The track with the most points in the area belongs to the takeoff or landing area. For the track that takes off from the same runway, the distance between its last track point and each corridor entrance is calculated. For the track that lands on the same runway, the distance between its first track point and each corridor entrance is calculated. The corridor entrance to which the track is closest belongs to.

Images