Keywords: Distance, trajectory, homotopy, energy, hierarchical clustering, airspace, major flows

In this paper, the problem of aircraft trajectories representation and analysis is addressed. In many operational situations, there is a need to have a value expressing how trajectories are close to each other. Some measures have been previously defined, mainly for trajectory prediction applications, all of them being based on distance computations at given positions in space and time. The approach presented here is to consider the trajectory as a whole object belonging to a functional space and to perform all computations in this space. An efficient algorithm for computing mathematical distance between trajectories is then presented and applied to the major flows extraction in the French airspace.

Future Air Traffic Management relies, in part, on the use of decision support tools (DST) to provide improved service to the user community under increasing traffic demand. Furthermore, this improvement has to be validated by the mean of system performance metrics such as complexity, robustness, capacity. As aircraft fly 4D trajectories, there is a strong need to quantify the associated trajectory accuracy in order to validate aircraft models and trackers. Such validation is usually based on a comparison between the actual trajectory and a reference by the mean of a trajectory distance. This last point is the key element of the whole process. Such trajectory distance is still needed for ATM applications and the goal of this paper is to present a new trajectory distance based on rigorous mathematical concepts. Although trajectories are well understood and studied, relatively little investigation on the precise comparison of trajectories is presented in the literature. A key issue in performance evaluation of ATM decision support tools (DST) is the distance metric that determines the similarity of trajectories. Most existing measures [START_REF] Jaynes | n open development environment for evaluation of video surveillance system[END_REF][START_REF] Senior | Appearance models for occlusion handling[END_REF] compute a mean distance of the corresponding posi-tions of two equal duration trajectories. Supplementary statistics such as variance, median, minimum, and maximum distances are also suggested to extend the description of similarity. In [START_REF] Needham | Performance evaluation metrics and statistics for positional tracker evaluation[END_REF], Needman proposed an alignment based distance metric that reveals the spatial transition and temporal shift between the given trajectories, and introduced an area based metric that calculates the total enclosed area between trajectories using trajectory intersection.

One main disadvantage of the existing approaches is that they are all limited to the equal duration (lifetime) trajectories. By duration we refer to the number of coordinate points that constitute the trajectory. These coordinates are sampled at different instances. Since the existing measures depend on the mutual coordinate correspondences, they cannot be applied to trajectories that have different durations. Conventional distance measures assume that the temporal sampling rates of trajectories are equal. They do not cope with the uneven sampling instances, i.e. varying temporal distance between the coordinates. Therefore, there is a need to develop other alternatives that can effectively measure the difference between unrestricted trajectories.There are a lot of ATM applications where such distance between trajectories is needed.

Aircraft model Inference All aircraft models are based on ODEs (Ordinary Differential Equation), including tabular ones (see Fig. 1). The aircraft model inference consists in answering the following question :Given a parametrized model and a goal trajectory, can we infer the best parameter values? A model can be viewed as a mapping from the control space into the trajectory space. The way to answer the previous question is then given by the closest model to the goal trajectory (see Fig. 2). In order to find the closest model in this trajectory γ Figure 2 Finding the best model from a given class. The green "grid" represents represents such class produced by the model and γ is the goal aircraft trajectory.

space, a reliable trajectory distance is needed. The model inference problem has to solve the accuracysmoothness dilemma :Over-fitted models are generally poor predictors. The previous construction gives the shortest path (and thus the distance) between the goal trajectory and the trajectory set which can be synthesized by the model.

Trajectory prediction

Air traffic management research and development has provided a substantial collection of decision support tools that provide automated conflict detection and resolution [START_REF] Erzberger | Conflict detection in the presence of prediction error[END_REF][START_REF] Brudnicki | User request evaluation tool (uret) conflict probe performance and benefits assessment[END_REF][START_REF] Vink | Eatchip medium term conflict detection: Part 1 eatchip context[END_REF],trial planning [START_REF] Mcnally | Field test evaluation of the ctas conflict prediction and trial planning capability[END_REF], controller advisories for metering and sequencing [START_REF] Swensen | Design and operational evaluation of the traffic management advisor at the forth worth air route traffic control center[END_REF][START_REF] Coppenbarger | Design and development of the enroute descent advisor (eda) for conflict-free arrival metering[END_REF], traf-fic load forecasting [START_REF] Meckiff | The tactical load smoother for multi-sector planning[END_REF][START_REF] Masalonis | Using probabilistic demand prediction for traffic flow management decision support[END_REF], weather impact assessment [START_REF] Kirk | Problem analysis resolution and ranking (parr) development and assessment[END_REF][START_REF] Sud | Air traffic flow management collaborative routing coordination tools[END_REF][START_REF] Evans | Reducing severe weather delays in congested airspace with weather support for tactical air traffic management[END_REF]. The ability to properly forecast future aircraft trajectories is central in many of those decision support tools. As a result, trajectory prediction (TP) and the treatment of trajectory prediction uncertainty is an active areas of research and development (eg [START_REF] Slattery | Trajectory synthesis for air traffic automation[END_REF][START_REF] Swierstra | Common trajectory prediction capability for decision support tools[END_REF][START_REF] Mondoloni | rajectory modeling accuracy for air traffic management decision support tools[END_REF][START_REF] Ryan | Review of trajectory accuracy methodology and comparison of error measure metrics[END_REF][START_REF] Mondoloni | Impact of factors, conditions and metrics on trajectory prediction accuracy[END_REF]).

Accuracy of TP is generally defined as point spatial accuracy (goal attainment) or as trajectory following accuracy. The last one can be rigorously defined by the mean of trajectory space. The first one is a limit case of the second by adding a weight function in the energy functional.

When we refer to trajectory prediction errors for a specific DST, we are typically comparing the predicted trajectory for a specific DST to the actual trajectory to be experienced by an aircraft. Discrepancies between these two types of trajectories typically affect the performance of the DST.

Radar tracker evaluation

The goal of a radar tracker is to eliminate the residual noise coming from the radars. It is a key element of the ATM system and its accuracy is one of the factors which determines the separation norm. In order to validate such trackers, an exact reference trajectory is generated and perturbed by a white Gaussian noise. This perturbed trajectory is then used as input of the tested tracker. The tracker generates an estimated trajectory which is compared to the reference trajectory. In order to do such comparison, a reliable trajectory distance is also needed.

Alternative route synthesis

Airspace congestion is related to aircraft located in the same area during the same period of time. Then, when congestion has to be minimized, algorithms have to separate aircraft in time (slot allocation), in space (route allocation) or both (bi-allocation). When route allocation is investigated, associated algorithms need alternative routes set in order to spread the traffic on them. A route is said to be alternative to another if it is different enough based on a trajectory distance.

Major flows definition When radar tracks are observed on a radar screen over a long period of time in a dense area, it is very easy to see major flows connecting major airports. The expression "major flows" is often used but never rigorously defined. Based on an exact trajectory distance and a learning classifier, it is possible to answer the following questions :Given a set of observed trajectories, can we spit it into "similar" trajectory classes? If yes, classes with highest number of elements will rigorously define the major flows. Given those classes and a new trajectory, can we tell if it belongs to a major flow and which one? The principle of the major flows definition is to use shape space to represent trajectory shapes as points and to use a shape distance. (the shape of a trajectory is the path followed by an aircraft, that is the projection in the 3D space of its 4D trajectory. The speed on the path has no impact). Major flows have not to be confused with highest density in the airspace. As a matter of fact, some approaches for major flows extraction consider the accumulated traffic in the airspace and build a kind of density map for which the highest areas are considered as major flows. This approach may be completely false as shown on figure 1. On this figure, four artificial trajectories share a common highest density area but, as it can be seen on the figure, no aircraft is flying this "high density trajectory".

Another approach consists in extracting major flows on a set of trajectories thank's to an efficient HMI and a bundling algorithm [START_REF] Marzuoli | Data visualization techniques for airspace flow modeling[END_REF]). The results produced by this kind of algorithm are quite similar to the ones presented in this paper but it is done manually and we propose to do it automatically.

As it has been shown in this section, mathematical distance between trajectories is a real need for many ATM applications. The next section of this paper presents some current trajectory distance metrics and shows their limitations. The third part gives a detailed mathematical description of our new trajectory distance. The fourth part introduces the associated algorithms implementation. Finally, the fifth part, presents the application of such algorithms to the major flows extraction of the French airspace with 8000 trajectories.

MATHEMATICAL DISTANCE BETWEEN TRA-JECTORIES

Introduction

In a vector space, distances are very well defined. If we consider two points P 1 = (x 1 , y 1 ) T and P 2 = (x 2 , y 2 ) T in a plane (see Fig. 4), the distance between them can be computed with the classical formula of the euclidean distance (see Fig. 4) : Figure 4 On the left, two points P 1 and P 2 has been drawn for which the classical Euclidean distance is shown in red. On the right, two trajectories are drawn (γ 1 , γ 2 ) for which one want to determine a mathematical distance.

d( P 1 , P 2 ) = (x 2 -x 1 ) 2 + (y 2 -y 1 ) 2 (1)
What is the distance, if now the points P 1 and P 2 are replaced by two trajectories γ 1 and γ 2 ? Trajectories are infinite dimension mathematical objects which are not easy to manipulate. We are looking for a mathematical distance between trajectories (γ 1 and γ 2 ) with the following properties :

• d(γ 1 (t), γ 2 (t)) = 0 ⇒ γ 1 (t) = γ 2 (t) • d(γ 1 (t), γ 2 (t)) = d(γ 2 (t), γ 1 (t)) • d(γ 1 (t), γ 2 (t)) + d(γ 2 (t), γ 3 (t)) ≥ d(γ 1 (t), γ 3 (t))
One of the main results of this paper is the establishment of such mathematical distance between aircraft trajectories.

Current Trajectory Distances

An aircraft trajectory is a time sequence of coordinates representing the aircraft path over a period of time and may be represented by a N-uple :T = {(x 1 , y 1 , z 1 , t 1 ), (x 2 , y 2 , z 2 , t 2 ), ..., (x N , y N , z N , t N )} where N is the duration.

The simplest metric used for computing the distance between a pair of trajectories is the mean of coordinate distance, which is given as

m 1 (T a , T b ) = 1 N N n=1 d n (2)
where the displacement between the positions is calculated using the Cartesian distance

d n = [(x a n -x b n ) 2 + (y a n -y b n ) 2 + (z a n -z b n ) 2 ] 1 2 (3)
Note that, the mean of distance metric makes three critical assumptions :

1. the durations of both trajectories are the same :

N a = N b = N 2. the coordinates are synchronized t a n = t b n 3. the time sampling rate is constant t a n+1 -t a n = t a m+1 -t a m
It is evident that the mean of distance is very sensitive to the partial mismatches and cannot deal with the distortions in time.

To provide more descriptive information, the second order statistics such as median, variance, minimum and maximum distance may be incorporated. For instance variance trajectory distance is defined as

m 2 (T a , T b ) = 1 N N n=1 (d n -m 1 (T a , T b )) 2 (4)
Although these statistics supply extra information, they inherit (even amplify) the shortcomings of the ordinary mean of distance metric m 1 . Besides, none of the above metrics is sufficient enough by itself to make an accurate assessment of the similarity.

Another possible candidate for the distance between two trajectories γ 1 and γ 2 will simply be to take the supremum norm (see Fig 5), that is : Since γ 1 and γ 2 are constant outside bounded intervals of R, the supremum is well defined. However, this metric is not sensitive to global properties of curves. In the Fig. 6, the curves γ 1 and γ 2 are at the same distance from γ 3 but have very different shapes. From an operational point of view, γ 1 is just a shifted For trajectories γ 1 , γ 2 with the same origin-destination pairs, γ 1 -γ 2 can be defined as a compactly supported mapping and an area distance between trajectories can be defined :

d ∞ (γ 1 , γ 2 ) = sup s∈R γ 1 (s) -γ 2 (s) (5)
d 2 (γ 1 , γ 2 ) = R γ 1 (t) -γ 2 (t) 2 dt 1 2 (6) d 2 ~area γ 1 (t) γ (t) 2 O D

Figure 7 Area distance between trajectories with the same origin-destination pairs

An extension of such area based distance metric is proposed in [START_REF] Needham | Performance evaluation metrics and statistics for positional tracker evaluation[END_REF]. The crossing points of two paths (where T a (p i ) = T b (p j )) are used to define regions Q j , j = 1, .., J between trajectories (see Fig. 8). For 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 111111111111111111111111111111111111111111111111 Figure 8 Area distance between trajectories with or without crossings.

1 γ 2 γ 1 γ 2 γ Q 1 Q 1 Q 2 Q 3
each region, a polygon model is generated and the enclosed area is found by the parameterized shape. The resulting distance is given by :

m 3 (T a , T b ) = N n=1 area(Q j ) (7)
This metric can handle more complex trajectories, however it is sensitive to entanglements of the trajectory, it discards the time continuity, and fails to distinguish two trajectories in opposite directions. Furthermore, it is not adapted to 3D trajectories. In order to introduce our new mathematical distance between trajectories, one must first give some representation definitions.

Representation

Since objects of interest are aircraft trajectories, we need to find an adapted framework in which computations may be made on trajectories as a whole. There are basically two ways of understanding what a trajectory is :

• The time/position approach. In this case, a trajectory can be represented as a mapping from a bounded interval of R (the life time of the trajectory) to R 3 or R 6 depending on whether speed is part of the data or not. Since there is an explicit dependence on time, there is a need to calibrate trajectories with time shifts for all applications involving trajectory comparison. We will see in the following that there is nevertheless a mean of reducing the problem so that origin of time is automatically calibrated.

• The shape approach. Here, trajectories are understood as paths and time is not directly relevant (from a more formal point of view, we take the quotient of the trajectories understood as mappings by the group of diffeomorphisms acting on time), so that we may assume that the underlying life time of trajectories is always the interval [0, 1]. This is the right framework for dealing with major flows estimation.

Trajectories as mappings

We will assume in the following that trajectories are given as mappings from a compact interval of R to R 3 . The case of mappings from R to R 6 (that is with explicit speed, for example as given by radar tracking filter) can be derived with minor changes and thus will not be addressed here. Since physical trajectories are smooth unless there is a perturbing noise, we made the choice to take all trajectories as smooth mappings from a compact interval of R to R 3 .

The first point to deal with is the necessary calibration of the origin of time for trajectories comparison. Remembering that there is an explicit dependence on time, one cannot just time shift one trajectory in time in order to make it coincident with another in order to compare them : this will result in forgetting distortions in time, that is trajectories with the same range (as mappings) but different positions at different times may become equal.

Since we choose to compare trajectories as mappings, a good candidate for computing the distance will be to integrate over time (like for the area distance) and to evaluate a mean error instead of the raw sum of squares :

d T (γ 1 , γ 2 ) 2 = 1 2T T -T γ 1 (t) -γ 2 (t) 2 dt (8)
with T > 0. Or, if we allow the mean to be weighted :

d T (γ 1 , γ 2 ) 2 = 1 T R h(t/T ) γ 1 (t) -γ 2 (t) 2 dt (9)
and h a positive summable function such that : R h(u)du = 1 [START_REF] Liao | Clustering of time series data -a survey[END_REF] This formula defines a semi-distance between trajectories γ 1 and γ 2 (see appendix A).

The previous family of semi-distances has nice features because of the scaling ability, but since it is not a single metric, it is difficult to use standard algorithms based on distances (for example, classification algorithms). There is thus a need for another definition of proximity between trajectories that will yield a single value while capturing interesting global characteristics.

Before introducing our homotopic distance between trajectories one must introduce how do we cope with time difference between trajectories.

Parametrization invariance

A very important constraint to take into account is the parametrization invariance: the shape of an object is independent on the way its contour is followed. In its seminal paper, Kendall introduced the notion of shape manifold [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF]: the originality of its work was the use of a differential geometry setting to implicitly enforce the invariance with respect to shape-preserving transformations. Curves were represented as finite sequences of distinguished points, called landmarks. Some related algorithms were eventually designed for air traffic analysis applications. In a study conducted by the Mitre corporation on behalf of the Federal Aviation Authority (FAA) [START_REF] Enriquez | Identifying temporally persistent flows in the terminal airspace via spectral clustering[END_REF], a spectral clustering algorithm was applied to sampled trajectories. Only the distance between landmarks was used, no invariance under euclidean transformations were imposed. Due to the high computational complexity, a random projection was first applied to the data in order to reduce the dimension of the samples. The most important limitation of this approach is that the shape of the trajectories is not taken into account when applying the clustering procedure unless a re-sampling procedure based on arc-length is applied: changing the time parametrization of the flight paths will induce a change in the classification. Methods based on times series as surveyed in [START_REF] Liao | Clustering of time series data -a survey[END_REF][START_REF] Rani | Recent techniques of clustering of time series data: A survey[END_REF] are appealing, but turn out to be inadequate for the present application. Finally, functional data statistics [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF][START_REF] Ramsay | Functional Data Analysis[END_REF] provides a powerful framework, still lacking the re-parametrization invariance. In this section, flight paths will be modeled as points in an infinite dimensional riemanian manifold. An intrinsic notion of distance exists in this setting and is defined as the infimum of the length of the paths connecting two points. Having this at hand allows the use of standard, distance based algorithms like k-means, k-mediods or hierarchical clustering.

Trajectories registration

A flight path may be modeled as a smooth curve γ : [a, b] → R 3 that maps a time to a position. Two distinct trajectories γ 1 , γ 2 are most of the time defined on different time intervals, say [a 1 , b 1 ] (resp. [a 2 , b 2 ]) for γ 1 (resp. γ 2 ), making the comparison between them quite awkward. This issue is well known in the field of functional data statistics as the registration problem. In a formal sense, it amounts to find a pair (φ 1 , φ 2 ) of strictly increasing diffeomorphisms

φ 1 : [0, 1] → [a 1 , b 1 ], φ 2 : [0, 1] → [a 2 , b 2 ] such that the transformed curves γ 1 • φ 1 , γ 2 • φ 2 ,
defined on the common interval [0, 1], are as similar as possible. The special problem instance:

min φ 1 ,φ 2 1 0 γ 1 • φ 1 (t) -γ 2 • φ 2 (t) 2 dt
gives the Fréchet distance between γ 1 , γ 2 . Computing the optimal φ 1 , φ 2 is a difficult task, unless the curves are assumed to be polygonal. Furthermore, as mentioned in [START_REF] Ramsay | Functional Data Analysis[END_REF], the registration procedure may remove some important features from the data: the extra degree of freedom provided by the so-called warping functions φ 1 , φ 2 may have the detrimental effect of registering curves that does not need it [START_REF] Ramsay | Curve registration[END_REF]. A discrete relative to the Fréchet distance is known as dynamic time warping and may be used to compare sampled sequences. Nevertheless, it suffers from the same drawback.

On the end of the other scale, a much simple procedure is to select only affine transformations for the warping functions. Given a trajectory γ : [a, b] → R 3 , the affine registration is γ • φ with:

φ : t ∈ [0, 1] → a + (b -a)t
It amounts to shift the time origin so as to make it coincident with 0, then to scale by the length ba of the time interval.

In between, registration procedures based on time landmarks or monotonic polynomial approximation may be used [START_REF] Ramsay | Curve registration[END_REF]. Most of the time, a penalty criterion must be added to the similarity measure in order to avoid the over-registration phenomenon. It worth mentioning a special procedure, that will be used in the sequel, that is more in line with geometry. Given a smooth curve γ : [a, b] → R 3 , its arclength is the smooth mapping:

s : t ∈ [a, b] → t a γ (u) du
The length l γ of the curve is just s(b). Assuming that γ never vanishes, s is strictly increasing, thus invertible. It induces a warping function:

ξ : t ∈ [0, 1] → s -1 (tl γ ) ∈ [a, b]
that is characterized by the property:

∀t ∈]0, 1[, D t γ • ξ(t) = l γ
where D t stands for the derivative with respect to t. This warping function is intimately related to the landmarks approach of [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF], as sampling evenly in the interval [0, 1] will result in a geometric even sampling on the curve itself (with respect to arclength). It will be denoted as the arclength warping in the sequel.

The manifold of paths

The idea of representing curves as point on an infinite dimensional manifold arises in the field of pattern recognition as an answer to the problem of assessing a degree of similarity between two shapes [START_REF] Peter | Riemannian geometries on spaces of plane curves[END_REF]. Within this frame, only closed curves were considered as they represent objects contours. In the context of air traffic, flight paths are never closed, unless the aircraft take off and land at the same airport, which is a quite uncommon for airliners. The initial mathematical model must be adapted to cope this specificity. For the sake of simplicity, all trajectories are assumed to be defined on the time interval [0, 1].

Definition 1

The space of immersions Imm([0, 1], R 3 ) is the set of smooth curves γ : [0, 1] → R 3 with nowhere vanishing derivative in the interval ]0, 1[. Generally speaking, an immersion will be a curve with nowhere vanishing derivative in the interior of its domain. It is clear that for such a curve the arclength is well defined and strictly increasing thus the geometric warping function exists. It may be used to perform a registration step to ensure that all curves are defined on [0, 1].

Given γ in Imm([0, 1], R 3 ), its derivative norm D t γ is a continuous mapping on the compact interval [0, 1] and thus has a non-zero minimum value m. If : [0, 1] → R 3 is a smooth mapping such that sup [0,1] D t < m, then γ+ will have a nowhere vanishing derivative and thus still belongs to Imm([0, 1], R 3 ). This indicates that this space has locally the structure of a vector space (in fact a Banach space) and globally the one of a differentiable manifold. To get rid of the influence of parametrization, the shape space is defined as a quotient with respect to all increasing diffeomorphisms of the interval [0, 1]:

E = Imm([0, 1], R 3 )/Diff + ([0, 1])
E inherits the manifold structure from Imm([0, 1], R 3 ). A point in E will be denoted by [γ] and is an equivalence class of mappings γ • φ with φ ∈ Diff + ([0, 1]). A tangent vector at [γ] is a couple ([γ], v) where v is a smooth mapping from [0, 1] to R 3 . This mapping must be understood as an infinitesimal displacement field on the base curve γ. As usual, the set of tangent vector is called the tangent bundle of E, denoted by T E. An riemanian metric can be introduced on E, in the spirit of [START_REF] Peter | Riemannian geometries on spaces of plane curves[END_REF]:

g [γ] (u, v) = 1 0 u(t), v(t) (1 + Aκ 2 (t)) D t γ(t) dt (11) + µ u(1), v(1) -µ u(0), v(0) (12)
where γ is a representative curve of [γ] and κ is the curvature of γ at t. The parameters A, µ are strictly positive real numbers that tune the respective importance of the curvature and the endpoints. The riemanian metric is invariant under a change of parametrization and thus does not depend on the particular choice of γ in the equivalence class [γ].

A smooth path between two points [γ 1 ], [γ 2 ] in E is represented by a smooth homotopy Φ ∈ Imm([0, 1], R 3 ), that is a smooth mapping from [0, 1] 2 to R 3 such that:

• Φ(0, •) = γ 1 (•), Φ(1, •) = γ 2 (•) • ∀(s, t) ∈ [0, 1] 2 , D t Φ(s, t) 0
The derivative of Φ with respect to the homotopy parameter s, denoted by D s Φn is a smooth curve on [0, 1], so that for a given s, the couple

([Φ(s, •)], D s Φ(s, •)) is a tangent vector in T E [phi(s,•)
] . An visual representation of an smooth homotopy along with the associated tangent vectors is given in Fig. 9.

(a, ) Using the riemanian metric [START_REF] Marzuoli | Data visualization techniques for airspace flow modeling[END_REF] on T E, the energy of a path Φ can be defined in the usual way:

E(Φ) = 1 0 g [φ(s,•)] (D s Φ(s, •), D s Φ(s, •)) ds (13)
It is equivalent for a path to minimize the energy or the length, the former is preferred as it saves a square root in the expression. The critical points of E are called geodesic paths. Since it is only a local condition (vanishing derivative), it may not correspond to a minimum of E. If such a global minimum exists, a path realizing it, is called a minimizing geodesic. In the finite dimensional setting, the Hopf-Rinov theorem may be invoked to prove the existence of a minimizing geodesic between arbitrary points. Unfortunately, it doesn't hold generally for infinite dimensional manifolds. It turns out that in the framework defined above, a minimizing geodesic exits between any two curves, thus making possible the definition of a distance on E:

For any couple ([γ 1 ], [γ 2 ]) in E 2 , the distance between [γ 1 ] and [γ 2 ] is given by:

d = 1 0 g [φ(s,•)] (D s Φ(s, •), D s Φ(s, •))ds (14)
where Φ is any homotopy between γ 1 , γ 2 realizing the minimum of E.

The distance d turns E into a metric space and can be used in any distance-based clustering algorithm.

ALGORITHM

Distance Algorithm

In order to compute the distance between two trajectories (γ 1 ,γ 2 ), a time regularization is first applied to both trajectories. Then, an homotopy Φ between γ 1 ,γ 2 is built for which a discrete grid is built in order to minimize its associate energy.

Let a be the origin of the trajectory γ. We have : γ(t) = a + t 0 γ (s)ds ,so a couple (a, γ ) (∈ W) with γ compactly supported defines a trajectory.

An homotopy between (a, γ 1 ) and (b, γ 2 ) is a continuous mapping Φ : [0, 1] → W such that Φ(0) = (a, γ 1 ) , Φ(1) = (b, γ 2 ). Intuitively, an homotopy is a continuous deformation between two trajectories.c

The deformation energy between γ 1 and γ 2 is linked to the distance between those trajectories and can be computed with the energy of the homotopy between γ 1 and γ 2 :

E(Φ) = 1 0 ∂v u ∂u 2 + R ∂γ u (s) ∂u 2 ds du (15)
In the case of a linear homotopy (which is the simplest one), the associated energy is given by :

Φ 0 (u, s) = [(1 -u).a + u.b] , (1 -u).γ 1 (s) + u.γ 2 (s) (16)
E(Φ 0 ) = b -a 2 + R γ 1 (s) -γ 2 (s) 2 ds (17)
There is an infinite number of homotopies shifting from γ 1 to γ 2 and our problem is to find the one with the minimum energy.

The deformation energy of a shape homotopy is obtained with a slight change in the expression for trajectories.

Figure 11 On this metric space each trajectory is represented by a point (blue point)..

E(Φ) = 1 0 ∂v u ∂u 2 + R ∂γ u (s) ∂u 2 . γ u (s) ds du (18)
In order to compute such energy, a grid on the homotopy connecting γ 1 to γ 2 is built, as shown on Fig. 10. This grid help us to compute an approximation of summation used in E(Φ). The optimization algorithm is searching for the z coordinate of each grid point in order to minimize E(Φ). One can show that such problem is convex (from the optimization theory point of view) and gradient like method can be used to find the associated minimum (quadratic programming has been used to solved this problem efficiently).

Clustering Algorithm

We consider a set of trajectories extracted from the radar track database of a given airspace. Having defined a distance between trajectories, one can gather together such trajectories in order to create clusters by using an adaptive clustering algorithm (hierarchical clustering). Such a clustering algorithm aims to partition the trajectory set into K clusters. To reach this goal, trajectories are consider as points in the associated metric space (see Fig. 11).

This algorithm uses two parameters, d min and d max , to respectively fuse clusters and create new clusters. Initially, each trajectory is considered as the centroid of a cluster. We then apply the three following principles one after the other:

• if two centroids are at a distance lower than d min , we fuse them into a single cluster, of which the resulting centroid is the barycenter of the two initial centroids. The barycenter is computed the following way :

µ i = 1 N i=N i=1 γ i (19)
• a new individual is aggregated to a cluster if its distance from the closest centroid is lower than d max and in this case we compute the new global centroid.

• Otherwise, we create a new cluster containing the single trajectory.

The number of clusters is also a result of the algorithm. An example of clustering resust is given on Fig. 12.

For each cluster c, one can compute also the following features :

• Number of trajectories in the cluster N c ;

• Mean trajectory which is the cluster centroid (γ c );

• Dispersion of the cluster:

N c i=1 γ j -γ c 2 ( 20
)
where . is the norm in the trajectory metric space.

The overall processing can be summarized by the Fig.

RESULTS

The algorithm has been applied to the French airspace with an heavy traffic of 8764 flights corresponding to June 27, 2015. This traffic has been extracted from the radar track database. Each trajectory being sampled every ten seconds, one has to manipulate about 5 million points, each on them having four coordinates (x, y, z, t). The traffic is represented on Fig. 14.

The initial step consists in computing the trajectory registration in order to remove the absolute time dependency. Then, the d min and d max distance have been fixed in order to apply the hierarchical clustering algorithm. Those distances have been establish during experimentations (d min = 2.30 d max = 4.5). Based on those distances, the hierarchical clustering algorithm has extracted 47 major flows for this day as shown on Fig. 15.

The algorithm has been implemented into C++ and executed on a IntelXeon3.2Ghz PC computer with an executing time of 30 seconds for extraction the major flow associated to the 8764 flights of June 27, 2015.

CONCLUSION

This paper has shown that distance between trajectories is a real need for ATM applications. Several ways of computing distances on the space of trajectories have been presented with their limitations. This family of metrics, scale based, is mainly useful for descriptive purpose and to quickly analyze a set of trajectories (for example, as a tool complementary to standard descriptive statistics).

Then, a concept originating from functional analysis has been introduced in order to work directly on trajectories as a whole. For more in depth analysis of trajectories, a new kind of distance has been introduced that is based on the energies of homotopies joining pairs of trajectories. This yield to a variational problem that cannot be solved directly, but may be reduced to a quadratic optimization problem. This kind of distance allows computations to be done on trajectories understood as shapes (or embeddings).

Based on this new distance between trajectories, an efficient major flows extraction has been developed with nice results on the French airspace for several thousands of trajectories. with T > 0. Or, if we allow the mean to be weighted :

d T (γ 1 , γ 2 ) 2 = 1 T R h(t/T )
((., .) denotes the usual scalar product) we have :

γ 1 (t) -γ 2 (t) 2 = ∆ -+ t -∞
2(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds with ∆ -= lim t→-∞ γ 1 (t) -γ 2 (t) 2 (this limit always exists since we have assumed at the beginning that trajectories are mappings from compact intervals). The expression of the mean square error becomes :

d T (γ 1 , γ 2 ) 2 = ∆ -+ 1 T R h(t/T ) t -∞
(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds dt By fubini theorem, we have :

d T (γ 1 , γ 2 ) 2 = ∆ -+ 1 T -T -∞
(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds+ (γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds This is in fact γ 1 (0) -γ 2 (0) 2 . γ 1 , γ 2 are constant mappings outside a compact interval of R, so both γ 1 and γ 2 are compactly supported. Furthermore the support of the mapping s → (γ 1 (s)-γ 2 (s), γ 1 (s)-γ 2 (s)) is included in the union of the supports of γ So the limit case T → +∞ is a convex combination of the initial and final differences (this will be 0 if trajectories are on the same origin destination, thus the limit of the d T is not a distance). Letting the weighting function h depend on time shift τ yields the final definition of a family of metrics :

d T,τ (γ 1 , γ 2 ) = sup τ∈R 1 T R h((t -τ)/T ) γ 1 (t) -γ 2 (t) 2 dt
(31) with the property that the limit case T → 0, T > 0 reduces to the supremum distance. This can be seen as a scale base distance, with T parameter being the scaling factor.

[

  EN-A-018] Trajectory Mathematical Distance Applied to Airspace Major Flows Extraction + + D. Delahaye * S. Puechmorel * S. Alam * * E. Feron * * * 1 Introduction

Figure 1

 1 Figure 1 General simulation model. The solver generate trajectories as close as posible to the reference trajectory thanks to the feedback path.

Figure 5

 5 Figure 5 Supremum norm distance

Figure 6

 6 Figure 6 Different trajectories with same sup distance

dvu vu initial pointFigure 9

 9 Figure 9 Smooth path between two curves

Figure 10

 10 Figure 10 Structure of the grid used for homotopy energy minimization. Each red point has 2D coordinates (x,y) for which an optimization algorithm is used for searching the z coordinates which minimize the energy of the homotopy connecting γ 1 to γ 2 .

13

 13

Figure 13

 13 Figure 13 Overall structure of the algorithm

Figure 15

 15 Figure 15 Major flows extracted from the hierarchical algorithm.

  )(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))dt ds and finally :d(γ 1 , γ 2 ) 2 = ∆ -+ s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds+ T -T (γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))( 1 s/T h(u)du)dsLimits as T → 0, T > 0 and T → +∞ can be obtained. If T → 0, T > 0, since 1 s/T h(u)du is bounded, we have :lim T →0 + d T (γ 1 , γ 2 ) 2 = ∆ -+ 0 -∞

= 1 0

 1 [START_REF] Brudnicki | User request evaluation tool (uret) conflict probe performance and benefits assessment[END_REF] and γ 2 , thus A > 0 exists such that this mapping vanishes outside [-A, +A]. This means that for T > A the first integral in the previous expression vanishes. The second has value : du)(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s)) ds[START_REF] Sud | Air traffic flow management collaborative routing coordination tools[END_REF] and here again, as soon as T > A this reduced to :du)(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s)) ds )du)(γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s)) dsby dominated convergence theorem, this gives :lim T →+∞ d T (γ 1 , γ 2 ) 2 = ∆ -+K A -A(γ 1 (s)-γ 2 (s), γ 1 (s)-γ 2 (s))ds[START_REF] Swierstra | Common trajectory prediction capability for decision support tools[END_REF] with K h(u)du. and finally, using again the fact that γ 1 , γ 2 are compactly supported, this is :lim T →+∞ d T (γ 1 , γ 2 ) 2 = ∆ -+K R (γ 1 (s)-γ 2 (s), γ 1 (s)-γ 2 (s))ds (28)Now, using the fact that : R (γ 1 (s) -γ 2 (s), γ 1 (s) -γ 2 (s))ds = ∆ + -∆ -(29) with ∆ + = lim s→+∞ γ 1 (s) -γ 2 (s) 2 , we have : lim T →+∞ d T (γ 1 , γ 2 ) 2 = (1 -K)∆ -+ K∆ + (30)

  γ 1 (t) -γ 2 (t) 2 dt (22)and h a positive summable function such that : Now, since the derivative of γ 1 (t) -γ 2 (t) 2 is : 2(γ 1 (t) -γ 2 (t), γ 1 (t) -γ 2 (t))

R

h(u)du = 1

[START_REF] Slattery | Trajectory synthesis for air traffic automation[END_REF]