Daniel Delahaye
email: delahaye@recherche.enac.frpuechmor@recherche.enac.fr

Stephane Puechmorel

Dsna

Air Traffic Complexity Map based on Non Linear Dynamical Systems

Keywords: Complexity, Dynamic Systems, Topological Entropy, Control Workload, Traffic Disorder

This work presents a new air traffic complexity metric based on non-linear dynamical systems. Previous work has shown that the structure and organization of traffic are important factors in the perception of the complexity of an air traffic situation. The new metric captures these important aspects of complexity by identifying the organization of trajectories in a traffic pattern. This paper investigates only the features of this new metric without quantifying directly the connection between complexity and the metric. Authors of previous work in this area have proposed metrics that generally have not explicitly addressed the effects of organization in the traffic flow on complexity. In order to capture the effect of organization, the metric is based on a dynamical system which fits as closely as possible the observations given by the aircraft positions and speeds. This approach, uses a non-linear dynamical system model that fits the observations without error. Based on this non linear dynamical system the method consist in the computation of the map of the associated Lyapunov exponents. This metric can be used to identify high (or low) complexity areas on a map, and, by capturing the organization properties of the traffic, captures some of the key factors involved in ATC complexity. Such maps are an example of the usefulness of these methods for comparing the relative complexity of different regions of airspace 1 .

Introduction

Air traffic control organizes air flows in order to ensure flight safety and to increase the capacity of the route network. Currently, about 7500 flights are registered everyday over France, which is a crossroad for the whole European airspace. This traffic generates a huge amount of control workload and the airspace is then divided into elementary sectors which are managed by air navigation controllers. For several years, a constant increase of air traffic has induced more and more congestion in the control sectors. Two strategies can then be applied to reduce such a congestion. The first one consists in adapting the demand to the existing capacity (slot-route allocation, collaborative decision making etc...). The second one adapts the capacity to the demand (modification of the air network, new design of the sectorisation, new airports, etc ...). For the two preceding approaches, the capacity of a sector is measured by the number of aircraft flying across the sector during a given period of time. The observation of controlled sectors shows that sometimes, controllers accept planes beyond the capacity threshold, while in other situations, they refuse traffic before the maximum capacity is reached. This phenomenon clearly shows that the operational metric alone cannot account for controller's workload. The goal of this study is to synthesize a traffic complexity indicator in order to better quantify the congestion in air sector, which will be more relevant than a simple number of aircraft which is independent of the traffic configuration. More precisely, our objective is to build some metrics of the intrinsic complexity of the traffic distribution in the airspace. Those metrics must capture the level of disorder (or organization) of any traffic distribution. Usually, metrics are focused on the speed vector distribution and the associated disorder metric captures only some features of the traffic complexity. The real objective of our work is to build a metric which measures the disorder or the organization of a set of trajectories in a 4D space (3D for the space and 1D for the time).

Those metrics are relevant for many applications in the air traffic management area. For instance, when a sectoring is designed [START_REF]Genetic algotithms for partitioning airspace[END_REF], the sectors have to be balanced from the congestion point of view and nowadays, only the number of aircraft is used to reach this objective. Another example where a congestion metric is needed is the traffic assignment [START_REF] Delahaye | Genetic algotithms for air traffic assignment[END_REF][START_REF] Delahaye | Aispace congestion smoothing by stochastic optimization[END_REF] for which an optimal time of departure and a route are searched for each aircraft in order to reduce the congestion in the airspace. Complexity metric may also be used to design new air networks, for dynamic sectoring concept, to define future ATM concepts (Free Flight) etc....Complexity metrics enable to qualify and quantify the performance of the Air Traffic service providers and enable a more objective consultation between airlines and providers.

The work presented in this paper is based on dynamical systems modeling of the air traffic. A dynamical system describes the evolution of a given state vector. If such a vector is given by the position of aircraft X = [x, y, z] T , a dynamical system associates a speed vector ˙ X = [v x , v y , v z ] T to each point in the airspace. The key idea is to find a dynamical system which modelizes the observed aircraft trajectories. Based on this dynamical system modeling, a trajectory disorder metric can be computed. In a first part, this paper will summarize the previous related works. The second part will present a linear dynamical system modeling for which the complexity metric can be represented into a complex coordinate system. In this system, it is very easy to identify any speed vector organization pattern. The third part introduces a non linear extension of the previous dynamical system modeling that may model any observed traffic. Such a non linear modeling can be used to produce maps of traffic complexity by identifying areas with high(low) complexity. This extension can also address time extension of the model and can then work straightly on the trajectory segments instead of speed vectors. The fourth part presents a tool developped by the DSNA which implements such algorithms for the computation of complexity map on real traffic situations.

Previous Related Works

The airspace complexity is related to both the structure of the traffic and the geometry of the airspace. Different efforts are underway to measure the whole complexity of the airspace.

Significant research interest in the concept of ATC complexity was generated by the "Free Flight" operational concept. Integral to Free Flight was the notion of dynamic density. Conceptually, dynamic density is a measure of ATC complexity that would be used to define situations that were so complex that centralized control was required [START_REF]Final report of rtca task force 3: Free flight implementation[END_REF].

Wyndemere inc. [START_REF] Inc | An evaluation of air traffic control complexity[END_REF] proposed a measure of the perceived complexity of an air traffic situation.This measure is related with the cognitive workload of the controller with or without knowledge of the intents of the aircraft. The metric is human oriented and is then very subjective.

Laudeman et al. from NASA [START_REF] Laudeman | Dynamic density : an air traffic management metric[END_REF] have developed a metric called "Dynamic Density" which is more quantitative than the previous ones and is based on the flow characteristics of the airspace. The "Dynamic Density" is a weighted sum of the traffic density (number of aircraft), the number of heading changes (> 15 degrees), the number of speed changes (>0.02 Mach), the number of altitude changes (>750 ft), the number of aircraft with 3-D Euclidean distance between 0-25 nautical miles, the number of conflicts predicted in 25-40 nautical miles. The parameters of the sums have been adjusted by showing different situations of traffic to several controllers. B.Sridhar from NASA [START_REF] Sridhar | Airspace complexity and its application in air traffic management[END_REF], has developed a model to predict the evolution of such a metric in the near future. Efforts to define "Dynamic Density" have identified the importance of a wide range of potential complexity factors, including structural considerations.

The traffic itself is not enough to describe the complexity associated with an airspace. A few previous studies have attempted to include structural consideration in complexity metrics, but have done so only to a restricted degree. For example, the Wyndemere Corporation proposed a metric that included a term based on the relationship between aircraft headings and dominant geometric axis in a sector [START_REF] Inc | An evaluation of air traffic control complexity[END_REF]. The importance of including structural consideration has been explicitly identified in recent work at Eurocontrol. In a study to identify complexity factors using judgment analysis, "Airspace Design" was identified as the second most important factor behind traffic volume [START_REF] Kirwan | Investigating complexity factors in u.k. air traffic management[END_REF].The impact of the structure on the controller workload can be found on the paper [START_REF] Histon | Introducing structural consideration into complexity metrics[END_REF][START_REF] Histon | Structural consideration and cognitive complexity in air traffic control, Proceedings of the 21 th Air Traffic Management for Commercial and Military Systems[END_REF]. Those papers show how strong the structure of the traffic (airways, sectors, etc...) is related with the control workload.

The previous models do not take into account the intrinsic traffic disorder which is related to the complexity. The first efforts related with disorder can be found in [7]. This paper introduces two classes of metrics which measure the disorder of a traffic pattern. The first class is based on geometrical properties and proposed new metrics which are able to extract features on the traffic complexity such as proximity (measures the level of aggregation of aircraft in the airspace), convergence (for close aircraft, this metric measures how strongly aircraft are closer to each other) and sensitivity (this metric measure how the relative distance between aircraft is sensible to the control maneuver). The second class is based on a dynamic system modeling of the air traffic and uses the topological entropy as a measure of disorder of the traffic pattern.

G.Aigoin has extended and refined the geometrical class by using a cluster based analysis [START_REF] Aigoin | Air traffic complexity modeling[END_REF]. Two aircraft are said to be in the same cluster if the product of their relative speed and their proximity (a function of the inverse of the relative distance) is above a threshold. For each cluster, a metric of relative dependence between aircraft is computed and the whole complexity of the cluster is then given by a weighted sum of the matrix norm. Those norms give an aggregated measure of the level of proximity of aircraft in clusters and the associated convergence. From the cluster matrix, it is also possible to compute the difficulty of a cluster (it measures how hard it is to solve a cluster). Multiple clusters can exist within a sector, and their interactions must also be taken into account. A measure of this interaction has been proposed by G.Aigoin [START_REF] Aigoin | Air traffic complexity modeling[END_REF]. This technique allows multiple metrics of complexity to be developed such as average complexity, maximum and minimum cluster complexities, and complexity speeds.

Another approach based on fractal dimension has been proposed by S.Mondoloni and D, Liang in [START_REF] Mondoloni | Airspace fractal dimension and applications[END_REF]. Fractal dimension is a metric comparing traffic configurations resulting from various operational concepts. It allows in particular to separate the complexity due to sectorization from the complexity due to traffic flow features. The dimension of geometrical figures is well-known: a line is of dimension 1, a rectangle of dimension 2, and so on. Fractal dimension is simply the extension of this concept to more complicated figures, whose dimension may not be an integer. The block count approach is a practical way of computing fractal dimensions: it consists in describing a given geometrical entity in a volume divided into blocks of linear dimension d and counting the number of blocks contained in the entity N .The fractal dimension D 0 of the entity is thus :

D 0 = lim d→0 log N log d (1)
The application of this concept to air route analysis consists in computing the fractal dimension of the geometrical figure composed of existing air routes. An analogy of air traffic with gas dynamics then shows a relation between fractal dimension and conflict rate (number of conflicts per hour for a given aircraft). Fractal dimension also provides information on the number of degrees of freedom used in the airspace: a higher fractal dimension indicates more degrees of freedom. This information is independent of sectorization and does not scale with traffic volume. Therefore, fractal dimension is a measure of the geometrical complexity of a traffic pattern. Some new geometrical metrics have been developed in [START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF]which are able to capture the level of disorder or the level of organization for some traffic patterns. For instance, in an artificial round about moving, the speed vector are very different even if the global moving is full organized without any changes in the relatives distance between aircraft. To capture those features, the covariance and the Koenig metrics have been developed. The first one is able to identify disorder or organization of translation movings. The second one identify organized curl moving.

All the previous metrics capture only one feature of the complexity and are not able to produce an aggregate metric which can capture all the possible situations (high-low density, how-low convergence, translation organization, round about organization etc ...). The topological entropy (Kolmogorov entropy) is the only metric which is able to capture most features of the complexity. The non-linear form is even able to identify any trajectories organizations (aircraft following the same path at the same speed). The next section will describe this metric in detail for the linear and the non-linear forms.

3 Linear Dynamical System Modeling

Principle

This metric is based on the modeling of a set of trajectories by a linear dynamical system. This enables to identify different structures of organization of the aircraft speed vectors such as translation, curve organizations, divergence, convergence or a mix of them. In the following of the paper, we will talk about a full organized traffic when the relative distances between aircraft do not change with time. For such a situation, the traffic is very predictable and very comfortable to address by a controller.

The key idea of this metric is to model the set of aircraft trajectories by a linear dynamical system which is defined by the following equation :

˙ X = A • X + B ( 2
)
where X is the state vector of the system :

X =   x y z   (3)
This equation associates a vector speed ˙ X to a position in the space coordinate X. This synthesis a particular vector field. The static behavior of this vector field is given by the vector B and the linear mapping between the speed vector ˙ X and the position vector X is given by the matrix A. The eigenvalues of the matrix A describe and summarize the evolution of the system. The real parts of those eigenvalues are related with the convergence or the divergence property of the system in the direction of the eigenvector. When such eigenvalues have positive real parts, the system is in expansion mode and when they are negative, the system is in contraction mode. Furthermore, the imaginary part of the eigenvalues are related with the level of curve tendency of the system. Depending of those eigenvalues, a dynamical system can evolve in contraction, expansion, rotation or a combination of those three modes.

The figure 1 give four typical examples for which the A matrix and the associated eigenvalues are computed. The small squares are the initial positions of aircraft at a given time(this represents the observation given by a radar for instance with the associated speed vector). As it can be seen on the figures, the aircraft are initially located on a circle with diameter 2d. The positions of aircraft 2 and 3 are symmetric to the x axis. The speed norms of aircraft are the same and have been fixed at v. In the four situations, the aircraft are evolving in a two dimensions space, so our state vector will be limited to those two dimensions :

X = x y (4)
This restriction enhance the simplicity of the equations for those examples, but the method works the same way in the three dimensional case. The computation of the four matrices A and vectors B is given by the following equations.

Situation 1 Situation 2 Situation 3

Situation 4

A = 0 0 0 0 A = v d -1 0 0 -1 A = v d +1 0 0 +1 A = v d 0 1 -1 0 B = v 0 B = 0 0 B = 0 0 B = 0 0 λ = v d 0 0 λ = v d -1 -1 λ = v d +1 +1 λ = v d +j -j (5)
It can be noticed that with such matrices, the equation2 fits exactly the observations for all situations. For instance, the position X of aircraft 1 in the fourth situation is given by :

X = -d 0 (6) So, ˙ X = A • X + B = v d 0 1 -1 0 . -d 0 + 0 0 = 0 v = v x v y (7)
which is the speed of the initial observation for this aircraft (3). For those four situations, the equation 2 fits exactly with the observations. This exact fitting is not always possible with a linear system but the LMS (Least Mean Square) procedure determines the system which fit better with the observation. The eigenvalues associated to the matrices A are given by the vectors λ and the corresponding loci are shown on the bottom part of figure 1. When the relative distances between aircraft remain unchanged with time (situation 1 and 4), it can be noticed that the real parts of the eigenvalues of matrix A are null. When the norms of the relative distances between aircraft diminish with time (situation 2) the real part of the eigenvalues are negatives; finally, when those relative distances increase with time (situation 3) the real part of the eigenvalue are positive. The largest those real part are in absolute value, the more those evolutions are fast.

The conclusion is then the following When the relative distances between aircraft change slowly with time (this mean that the relative speeds between aircraft are close to zero and the traffic has no interaction), the eigenvalues of matrix A of the associated dynamical system have a small real part (in absolute value).

In the complex coordinate system, one can identify a small vertical strip where are located the eigenvalues of the "organized" dynamical systems (see figure 2). It must be noticed that the LMS regression of the A matrix and the B vector is done with the help of a singular value decomposition for which the singular value are check in order to avoid numerical problem. As a matter of fact, when such a value are not large enough, the associated dimension will not be taken into account in the following computation, this will ensure that the matrix A can always be determinate for the significant dimension. 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 [START_REF] Inc | An evaluation of air traffic control complexity[END_REF] In this coordinate system, the organized situations (in translation, curve or both) are located around the imaginary axis (shadow area). which represent the eigenvalues with a small real part.

A dynamical system can also be considered as a map T from the state space X to itself. Based on this mapping T , one can define the topological entropy (or Kolmogorov entropy) of the dynamical system which measure the level of mixing of X by T . This entropy is associated with the changes of the relative distances between points from X by T . The topological entropy of a dynamical system is then a disorder indicator of the distribution of the aircraft in the considered airspace. This entropy is not a statistical metric and is more adapted to the air traffic control system for which there are few aircraft in the sectors (statistical entropy needs a lot of samples in order to be reliable).

Kolmogorov Entropy Computation

The Kolmogorov entropy is computed with the help of the eigenvalue decomposition of the matrix A :

A = P -1 • D • P ( 8
)
where D is the diagonal matrix of the eigenvalues. Based on the observations of the aircraft (positions and speed vectors), the dynamical system has to be adjusted with the minimum error. This fitting has been done with a Least Square Minimization (LMS) method [START_REF] Bjorck | Numerical methods for least squares problems[END_REF]. For each considered aircraft i, it is supposed that position X i = [x i , y i , z i ] T and speed vector

V i = [v x i , v y i , v z i ] T are
given. An error criterion between the dynamical system model and the observation is computed :

E = i=N i=1 V i -A • X i + B (9)
The minimization of this criterion with the LMS induces an optimum set of parameters : A = P -1 • D • P and B. The eigenvalues of the matrix A (matrix D) are complex numbers. In order to produce a scalar metric, the larger absolute value of the real part of those complex eigenvalues is computed. The sign of this larger real part is related with the mode of the system (contraction or expansion). If those eigenvalues are drawn in the complex coordinates, they give information about the level of contraction/expansion of the system and also its curve tendency. Such a metric can identify any speed vector organization (translation, curve or both) because it is sensitive to relative distances between aircraft. The analysis of such loci may be used in order to compare several traffic situation. Linear dynamical system modeling enables to produce an aggregate metric associated to any traffic situation and can recognize any global organization pattern. If detail metric is needed in order to identify high (low) complexity areas, a non-linear dynamical system has to be used.

Non Linear Extension 4.1 Introduction

Evolution equations of the form :

˙ X(t) = V ( X(t), ω) (10)
with V C 2 vector field depending on parameters ω, describe systems which integral curves may fit exactly the observed trajectories (which is not generally possible with linear models).

There are many classical ways of obtaining a class of parametrized vector fields which fulfill the fitting requirement. Among them, vector splines allow a control on the smoothness of vector fields, which is important in our case since civil aircraft maneuvers are based on low acceleration guidance laws. Like in the linear case, the approximating model will be first computed, then Kolmogorov entropy will be estimated using sampled integral curves of the underlying vector field.

Div-Curl Splines

Computing topological entropy for a given traffic situation requires interpolating a vector field given only samples (positions and speeds of aircraft at a given time). Vector spline interpolation seeks the minimum of a functional of the form :

1 2 D L V ( X) 2 d X + 1 2 m i=1 V ( X i ) -V i 2 ( 11
)
where V is a vector field defined on a domain D ⊂ R n , L is an elliptic differential operator controlling smoothness of the solution and ( X i , V i ) i=1...m are the interpolation data [START_REF] Suter | Fast evaluation of vector splines in three dimensions[END_REF]. By introducing the adjoint operator L T , optimal vector field can be shown to be a linear combination of shifted version of the elementary solution kernel of the differential operator L T L.

A special case is the so-called "div-curl" splines with the criterion :

R 2 α ∇div V ( X) 2 + β ∇curl V ( X) d X (12)
with α, β positive reals controlling the smoothness of the approximation by focusing on constant divergence or constant curl.

Dynamic Div-Curl Splines

Div-Curl splines work well for spatial vector field interpolation (wind or marine flows data for example). However, in our case, aircraft trajectories may intersect even if there is no conflict. In such situations, the classical div-curl interpolation cannot provide a solution since the crossing point will become a singular point of the approximating vector field. To deal with that we introduce a family of spatio-temporal splines obtained as optimal solution of a dynamic interpolation problem, namely find a time-dependent vector field V (t, X) defined on [0, T ] × D and continuously differentiable up to order 2 in time coordinate that minimize :

T 0 D ∂ V (t, X) ∂t 2 + α L V (t, X) 2 dxdt (13)
under the constraints :

V (t i , X i ) = V i , i = 1 . . . m (14)
where α is a positive real. This value controls the relative importance of the vector field variation through time over the discrepancy of V as measured by the differential operator L.

Taking α = 0 will yield to a constant vector field over time, while α → +∞ will focus on the differential part.

Solving the associated variational problem

Let L be an elliptic differential operator of order p with constant coefficients. L t will denote the adjoint operator of L. We will assume in the following that the vector field V has fixed value at 0 and T and that for all t ∈ [0, T ], V (t, .) belongs to the Sobolev space H p . Furthermore, the mapping t → V (t, .) must be at least two times continuously differentiable.

Let γ : [0, T ] → C p (R n , R n ) a variation of X (that is time-dependent H p vector field γ(0, .) = γ(T, .) = 0). Simple calculus shows that the gradient of :

T 0 D ∂ V (t, X) ∂t 2 + α L V (t, X) 2 d Xdt (15)
with respect to γ, is the linear mapping :

γ → T 0 D γ, ∂ 2 V ∂t 2 -αL T L V d Xdt (16)
The constraints on X impose that :

γ(t i , .) = 0, i = 1 . . . m (17)
One can find an elementary solution e of the partial differential equation

∂ 2 V ∂t 2 -αL T L V (18)
and use its translate e t i , X i to represent the evaluation operator at (t i , X i ). Combining this with the gradient previously obtained and using Lagrange multipliers (denoted λ i , i = 1 . . . m) shows that the optimal solution is of the form :

V = m i=1 λ i e t i , X i + V 0 (19)
with V 0 an element of the kernel of the operator :

∂ 2 V ∂t 2 -αL T L V (20)

Complexity Map Computation

The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of the underlying dynamical system called Lyapunov exponents. In order to figure out what Lyapunov exponents are, let consider a point and look at its evolution when transported by the dynamical system.

Let X be fixed (initial point) and let φ be a point trajectory of the dynamical system given by :

φ(t, X) = X + t 0 V φ(s, X)ds (21)
Assume now that trajectory is disturbed by a small perturbation , we have :

φ(t, X + ) = φ(t, X) + D X φ(t, X) • + o( ) (22)
where D X φ(t, X) is the differential of the vector field at X that satisfies :

D X φ(t, X) dt = D X V (φ(t, X)) • D X φ(t, X) (23)
The Lyapunov exponents are the singular values of the matrix D X φ(t, X) and can be thought as local shear values for the dynamical system. When Lyapunov exponents are high, the trajectory of a point under the action of the dynamical system is very sensitive to initial conditions (or parameters on which the vector field may depend), so that situation in the future is unpredictable. On the other hand, small values of the Lyapunov exponents mean that the future is highly predictable (expected to be comfortable for a controller). So, the Lyapunov exponent map determines the area where the underlying dynamical system is organized. It identifies the places where the relative distances between aircraft do not change with time (low real value) and the ones where such distance change a lot (hight real value).

Let us now describe the practical procedure for computing complexity maps. First of all, the optimal dynamic div-curl approximation for the observed trajectories is computed, based on the defining equations. That step requires a linear system solving.

The second step computes the second derivatives matrix D X φ at each point of the grid for φ trajectory starting at X. This is done by solving the differential equation 23 with a Runge-Kutta integrator. The complexity value at point X is then obtained by averaging Lyapunov exponents over the time :

κ( X) = 1 n i=n i=1 D X V (φ(t, X)) 2 (24)

Results

In order to validate the previous method, complexity maps have been computed for several senarii and have produced relevant results. Those maps have been computed by a DSNA software called CASTOR for Complexity for Airspace and Traffic ORganization. The input data of this software are trajectory samples coming from radar or any traffic simulator. Those trajectory samples are then used to compute several complexity metrics which are summarized by maps on an HMI. The structure of this project is given on figure 3. This application is network oriented and the main components are the following :

• Traffic database. Berkeley Data Base Manager(DBM) has been used for trajectory storage. The trajectory samples are summarized by the following field :id,time,lat,long, altitude, Vx,Vy,Vz. This database may be fed by COURAGE traffic simulator, radar simulators, real radar samples. The database implements a time-space selection for trajectory extraction and can control the data flow by the mean of a time accelerator system. This database has been implemented as a network traffic server (see figure 4).

• Complexity computation. Several complexity indicators may be selected in order to build the associated 3D map. The current version of the software proposes the following indicators : proximity, convergence, covariance, koenig transation, koenig rotation, cluster and topological entropy. This last indicator is the one which is based on the Lyapunov exponents computation (see previous section). This computation has been implemented as a computing server (see figure 5).

• Human Machine Interface. This HMI is given on figure 6 where a map of Europe is drawn. The traffic evolves on this map and the user may select traffic sample in 3D. The complexity indicator is selected on the upper left part of this HMI. The separation norm being different in the horizontal direction (ex 5 nmi) and the vertical direction (1000ft), the complexity has been computed with a non-isotrop distance, meaning that two aircraft separated by 5 nmi in the horizontal plane or 1000ft in the vertical direction have the same "distance" in the new non-isotropic coordinate system for which the unit are separation norm in any direction. A first test has been used in order to validate the computation of indicators. This first situation is given on figure 7. This first test gather together four artificial traffic situations each involving eight aircraft. On the bottom left the aircraft are full organized in parallel; on the upper left those eight aircraft are converging on the same point; on the upper right the aircraft follow a circle trajectory; and on the bottom right a random traffic situation has been created. A first complexity map has been computed on this traffic situation for which the proximity indicator has been selected (see figure 8). This indicator identify areas where traffic is more aggregated. This map shows that the four situations present the same proximity feature. A second map is presented on figure 9 for which the convergence indicator has been selected. This indicator shows areas where aircraft are converging.This map presents complexity for the random traffic (bottom right) and the convergent traffic (upper left). Finally, the figure 10 shows the associated Lyapunov exponents map for which only the full organized traffic (bottom left) do not generate complexity. It can be noticed that the converging situation (upper left) do not generate much more complexity than the traffic organized on a circle (upper right). This is absolutely normal in the sense the system knows that the converging situation may be easily solved by putting all the aircraft on the round about system (right or left). So, it is easy to shift from a full convergent situation to a situation where all the aircraft are well separated. The second test presents full organized miles in trail trajectories (from south west to north east) with two random traffic situations (see figure 11). All the classical geometrical indicators do not exhibit this miles in trail organization as it can be seen for instance on the figure 12 for which the proximity indicator has been selected. Only the Lyapunov exponent are able to identify such trajectory organization (see figure 13). This figure shows clearly a complexity valley on the miles in trail direction. This organization may have been detected even if the miles in trail trajectories would have been structured on a curve trajectory. That is the strong point of this metric: Lyapunov exponents are able to identify any kind of trajectory organization.

All those indicators may be computed on real traffic samples over France as it is shown on figure 14.

Conclusion

We have presented in this paper a new air traffic complexity metric based on non-linear vector field model of air traffic. Extending previous results on topological entropy, this method allows identification of traffic pattern organization in its full generality while previous works were limited to specific aspects of the complexity. Furthermore, since we are processing trajectory segments instead of samples at a given time, the induced metric is more adapted to the air traffic management which works mainly with trajectories. This facility enables to take the flight plan into account in the computation in order to figure out the complexity in the near future like controllers are doing when they address the traffic. The quasi-interpolation

Figure 1 :

 1 Figure 1: Eigenvalues loci for several typical situations

Figure 2 :

 2 Figure 2: Impact of the eigenvalues on the dynamic of a system

Figure 3 :

 3 Figure 3: CASTOR project structure

Figure 5 :

 5 Figure 4: Database server structure

Figure 6 :

 6 Figure 6: HMI structure

Figure 7 :

 7 Figure 7: First test

Figure 9 :Figure 10 :

 910 Figure 9: Convergence map

Figure 11 :

 11 Figure 11: Second test situation

Figure 12 :

 12 Figure 12: Proximity indicator on the second situation

Figure 14 :

 14 Figure 14: Example of map computation for a traffic sample over France

The complexity software described in this paper has been developped by NeoMetSys for DSNA. We may thank Ren Zanni for its advices and expertise on operational ATM system