David Gianazza

Kévin Guittet

SELECTION AND EVALUATION OF AIR TRAFFIC COMPLEXITY METRICS

This paper presents an original method to evaluate air traffic complexity metrics. In previous works, we applied a Principal Component Analysis (PCA) to find the correlations among a set of 27 complexity indicators found in the literature. Neural networks were then used to find a relationship between the components and the actual airspace sector configurations. Assuming that the decisions to group or split sectors are somewhat related to the controllers workload, this method allowed us to identify which components were significantly related to the actual workload. We now focus on the subset of complexity indicators issued from these components, and use neural networks to find a simple relationship between these indicators and the sector status.

Introduction

How complex is it to control a given air traffic situation? Over the last years, this question has become more and more crucial, as air traffic increased, thus creating bottlenecks in the air traffic control system. There arised a need for performance and complexity metrics, allowing to assess how the current system is operated and how it could be improved.

As a consequence, air traffic complexity has become the subject of many studies -often in relation to the controller's workload -and a multitude of complexity metrics have been proposed. The inadequacy of the aircraft count to appropriately reflect the traffic complexity has now been acknowledged for a long-time, and complementary indicators such as "traffic mix", "number of potential conflicts" and others, have been (and still are) designed. A linear combination of these variables, often referred to as dynamic density, is likely to better fit traffic complexity than individual indicators. It is used throughout most studies, where the correllation of a set of indicators with a quantifiable variable (subjective ratings, number of interactions with the computer,...) assumed to represent the actual traffic complexity, is maximized.

A possible shortcoming of this methodology is that potentially non-linear relations between indicators are missed (see [START_REF] Hilburn | Toward a non-linear approach to modeling air traffic complexity[END_REF] and the concern of Eurocontrol when writing calls for proposals). But, more importantly, the choice of the dependent variable is crucial to determine how well complexity is actually measured. Indeed, physical activity, as used in [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF] and [START_REF] Majumdar | The factors affecting airspace capacity in europe: A framework methodology based on cross sectional time-series analysis using simulated controller workload data[END_REF], miss the important cognitive part of the controller activity. On the other hand, physiological indicators ( [START_REF] Crump | Review of stress in air traffic control: Its measurement and effects[END_REF], [START_REF] Averty | Evaluating a new index of mental workload in real atc situation using psychological measures[END_REF]) seem difficult to exploit and how well they relate to traffic complexity is unclear. Finally, widely used subjective ratings ( [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]) provide high quality data (as they obviously relate to the kind of complexity investigated), but are often seen as subject to biases (such as the recency effect denounced in [START_REF] Averty | Evaluating a new index of mental workload in real atc situation using psychological measures[END_REF], and the possibility of raters errors in the case of "over-the shoulder workload ratings" [START_REF] Mannings | The relationship between air traffic control events and measures of controller taskload and workload[END_REF]). In all of these cases, data are very expensive to collect, as they require the active participation of controllers. Databases are often small and might exhibit low variability, which may in turn harm the statistical relevance of the results. This phenomenom is acknowledged in [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], where the overfitting of data is clearly a consequence of a lack of observations rather than a misspecification of the neural network. Finally, as these complexity metrics may be used to design computer-assisted control tools or traffic management tools, and to help organizing airspace, it is surprising to notice that the question of the relevance of the complexity measured to the final goal is scarcely discussed. The question really is to understand which complexity is measured and how well it relates to the foreseen application (benchmarking, improvements in airspace organization, design of new tools...).

Our research is initially motivated by former studies on optimal airspace sector configurations ( [START_REF] Gianazza | Optimization of air traffic control sector configurations using tree search methods and genetic algorithms[END_REF], [START_REF] Gianazza | Optimal combinations of air traffic control sectors using classical and stochastic methods[END_REF]) and intends to improve the criterion used therein to evaluate sector configurations. The basic idea, introduced in [START_REF] Gianazza | Optimisation des flux de trafic aérien[END_REF], is that the decisions to split a sector, mostly taken when the controller is close to overload, are linked to traffic complexity and may therefore provide an acceptable dependent variable. Interestingly enough, collecting data on sector configurations does not require controllers active participation, as current outcomes from control centers can be used, while related flight informations are available from recorded radar tracks. As such, raw data needed in our study are noticeably cheap to collect and might be produced in large quantities. The price to pay is that these data are noisy, as we may not be sure that a sector splitting (resp. merging) decision is directly related to overload (resp. underload). Other factors might distort data, such as training of unexperienced controllers, meteorological hazards, military airspace use... However, we will assume that the impact of these phenomena on the accuracy of the results is limited, particularly because of the kind of complexity we are looking at here. Indeed, this work is conducted in the perspective of future pre-tactical applications (e.g. sector planning) and thus does not ask for as much details as studies of instantaneous workload would (on the opposite, benchmarking of ATC centers would require an even coarser granularity, as indicators are averaged on wide temporal and geographical horizons [START_REF] Flynn | Traffic complexity indicators and sector typology analysis of u.s. and european centres[END_REF], [START_REF]Note de synthèse sur l'indicateur de complexité pru[END_REF]).

This paper is a continuation of [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF], which started with a principal component analysis on a set of indicators found in the literature, and where we proposed an original method allowing to select which components were actually related to the controller's workload, considering the status of the control sector. Neural networks were used to investigate the link between complexity indicators and sector configurations, as non-linear interactions were suspected.

For practical purposes, we would rather avoid to compute all the indicators and then apply a transformation matrix in order to find the relevant components. So, our aim now is to predict the sector status directly from a small subset of relevant complexity indicators, issued from the previous selection of principal components.

The paper is organized as follows. Section I briefly describes the indicators used throughout the study, while section II presents the raw data from which the final database is built. Section III summarizes the results of the Principal Component Analysis (PCA) which was performed in [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF]. Neural networks are introduced in section IV. The results on the complexity metrics selection and evaluation are presented and discussed in sectionV. Section VI concludes.

I. Air traffic complexity indicators

TABLE I CHOSEN SUBSET OF AIR TRAFFIC COMPLEXITY INDICATORS

The accuracy of the results of a study related to air traffic complexity is strongly dependent of the diversity and quality of the chosen individual complexity indicators. Many have been suggested to help describe the controllers workload, and it is hardly possible to implement the entire pool. In order to limit the number of variables to be (re)programmed and present indicators that are representative of the dynamic density literature, we focused on the ones selected by Kopardekar [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] in its unified complexity metric 1 . These indicators, and references to studies where they were used and where definitions may be found, are presented in Table I. We also implemented several indicators inspired by studies conducted elsewhere in the SDER (former CENA). Definitions are indicated in appendix A 2 . Finally, we also used incoming flows as explanatory variables, as they may be a significant factor in the decision to split (or merge) a sector.

II. Input data

The indicators are computed every round minute of the day, using recorded radar data, environnement data (sector description), and recorded sector configurations of the five french ATC centers. The sector configurations are recorded every round minute of the day, which explains our choice concerning the frequency at which we compute the indicators.

Radar data is available in several forms: records made by each center, with one position every twelve seconds, in average, and a global record of the five centers, with one position every three minutes. Several months of global records were available, whereas the centers local records were not readily available, at least for a sufficiently long period of time. So we used the global records (made by the IMAGE system), and interpolated the aircraft positions in order to get one position per minute. As many trajectory changes may occur within three minutes of flight, the computed positions are not highly accurate, and this may introduce a bias in the indicators values. However, this bias is most probably of small importance in our problem: we just want to predict when a sector will be merged into another one, or split in several smaller sectors. We are not considering the instant workload, which may require a very high level of accuracy on the aircraft position, speed, and so on. To be sure that this bias is small, we should compare the computed positions, and maybe also the indicators values, 1 Though we were not always able to find an explicit formula, and thus missed seamingly important indicators like, e.g., "MET_airspace structure". Note that this difficulty to get clear definitions is also reported by Eurocontrol in [START_REF]Cognitive complexity in air traffic control, a litterature review[END_REF]. 2 Further informations and discussions about indicators are to be found in the internal note [START_REF] Guittet | Analyse descriptive des indicateurs de complexité du trafic aérien à partir des données image et courage[END_REF].

using local centers records, and global records, on small data samples. This is left for future work.

Several months of recorded traffic are available. However, considering the volume of data, it would be tedious to run several experimentations on very large data samples. So, we have restricted our choice, at least for the moment, to one day of traffic (1st june, 2003). Once we have found the most significative complexity indicators, it will be possible to re-train the neural network on larger data samples.

On the chosen day, 103 different sectors were armed. The term "sector" means here either an elementary sector, or a set of elementary sectors merged together, and handled on a single controller's working position. The air traffic complexity indicators were computed for each of these sectors, every minute of the day, together with the sector status (merged, armed,o rsplit). This data was split into two sets : about sixty percent was randomly selected in order to train the neural network, and the rest was used to test the trained network on fresh data.

This single day of traffic already provides a big volume of data, as detailed in table II, with a great diversity of geographic sectors, and with enough data in each class of sector status.

III. Principal component analysis

Before applying the neural network to complexity indicators and sector statuses, let us first give a brief summary of the principal component analysis results. The reader may refer to [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF] for the details.

The aim of the principal component analysis is to find the colinearities between the indicators, and thus reduce the dimensionality of our problem. From the initial 27 complexity indicators, 6 main components were identified (corresponding to eigenvalues greater than 1), that covered more than 76 % of the variance of the data set.

These components are interpreted as follows: C 1 may be seen as a "size factor" and is strongly representative of the aircraft count, C 2 is related to the ground speed variance, and the aircraft vertical evolutions, C 3 is highly correlated with incoming flows, C 4 seems mostly related to converging flows and anticipation of conflicts, C 5 is linked with divergent flows, and C 6 is strongly correlated with the vertical proximity measures ( [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]), and could stand for the monitoring of vertical separation (near the minimas).

Notice that we extracted only 6 components, thus significatively less than the 12 components (briefly) described in [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF]. This might be explained by the lack, at the point of the project, of indicators related to the sector geometry.

Component Variables

C2 avg_vs,

σgs gs , σ 2 gs , N cl , N ds C3 F5, F15, F30, F60 C4 inter_hori, insen_c, Conv, creed_ok, creed_pb C5 insen_d, Div C6 vprox_1, vprox_2

TABLE III VARIABLES BELONGING TO THE PRINCIPAL COMPONENTS

Table III lists the indicators providing the best explanation for each component between C 2 and C 6 .F o rC 1 , we will simply use the number of aircraft Nb, as it is the most representative of the overall size of the problem. For C 2 and C 3 , the selected variables have an absolute value of correlation above 0.4 (that is the correlation with the component axis). For these two components, the chosen variables are highly representative of the corresponding components. The choice is less easy for component C 4 , where several indicators with the highest correlation values are, in fact, not very meaningful, as they are already involved in previous components. This means that we probably should refine again the definition of some of the indicators, in order to have a clearer interpretation of this component. For the time being, we have selected the variables with the highest absolute values of correlation, and which did not appear in previous components. For C 5 and C 6 , the selected variables are well correlated with the component (0.4 to 0.7).

IV. Neural networks

A. Presentation

Artificial neural networks are algorithms inspired from the biological neurons and synaptic links. An artificial neural network is a graph, with vertices (neurons, or units) and edges (connections) between vertices. There are many types of such networks, associated to a wide range of applications: pattern recognition (see [START_REF] Bishop | Neural networks for pattern recognition[END_REF] and [START_REF] Ripley | Pattern recognition and neural networks[END_REF]), control theory,... Beyond the similarities with the biological model, an artificial neural network may be viewed as a statistical processor, making probabilistic assumptions about data ( [START_REF] Jordan | Neural Networks[END_REF]). Some train data is used to determine a statistical model of the process which produced this data. Once correctly trained, the neural network uses this model to make predictions on new data. Neural networks are closely related to the Bayesian probabilities. They may be used for unsupervised learning (density estimation problems), and, mainly, supervised learning problems (regression, classification). Density estimation is not in the scope of our paper, so we will not detail it. The aim of regression is to find a statistical model producing an output y from input variables (let us denote them by x), so that the output y is as close as possible to a target variable, which we shall denote by t. In the case of classification problems, the target variables represent class labels, and the aim is to assign each input vector x to a class.

We will use a specific class of neural networks, referred to as feed-forward networks, or multilayer perceptrons (when the activation function is logistic). In such networks, the units (neurons) are arranged in fully-connected layers: an input layer, one or several hidden layers, and an output layer. Figure 1 shows an example of such a network.

For our problem, we have chosen three-layers feed-forward networks, denoted I α H β O γ in the rest of the paper, with α units in the input layer, β units in the hidden layer, and γ units in the output layer. The input variables are normalized, by substracting the mean value and dividing by the standard deviation.

The output vector y =( y 1 , ..., y k , ..., y q ) T is expressed as a function of the input vector x = (x 1 , ..., x i , ..., x p ) T as follows:

y k =Ψ( q j=1 w jk Φ( p i=1 w ij x i + w 0j )+w 0k ) (1)
where the w ij and w jk are weights assigned to the connections between the input layer and the hidden layer, and between the hidden layer and the output layer, respectively, and where w 0 j and w 0 k are biases (or threshold values in the activation of a unit). Φ is the sigmoid logistic function :

Φ(z)= 1 1+e -z (2)
For the transfer function Ψ applied to the output layer, we have chosen the softmax function:

Ψ(z k )= e zk C m=1 e zm (3)
With this choice for Ψ, the error function minimized during training is the following loglikelihood, known as cross entropy:

E(w)=- N n=1 C k=1 t (n) k ln(y (n) k ) (4)
where t (n) and y (n) are the n th target and output vectors, respectively.

The error is a function of the network's parameters w (weights and biases). The training aims at choosing these parameters, so as to minimize the chosen output error function. As we adress a classification problem -assign each input vector (a list of complexity indicators values) to a class representing the sector status (merged, armed,o rsplit) -it is best to minimize the cross-entropy instead of the quadratic error often used in regression.

The neural network is designed with one output unit per class (merged, armed,o rsplit). A target vector t (n) with value (1, 0, 0) means that the considered sector was merged with other sectors when the n th measure of the vector of complexity indicators was made. Armed sectors will be represented by (0, 1, 0). A value of (0, 0, 1) means that the sector was previously split in two or more sectors at the time x was measured. Of course, the actual output y (n) of the neural network will not be made of exact values 0 or 1. It will be a triple (a, b, c) of floatingpoint values between 0 and 1, each value being the probability to belong to the corresponding class. The input vector x (n) will be assigned to the class of highest probability.

The number of input units of the network is the cardinal of the evaluated set of metrics. For example, if we consider {V, Nb} we will use a network with two input units. The hidden layer comprises 15 units.

The nnet package of the R language was used (see http://www.r-project.org/ for details on the R language and environnment). In this package, developped by Pr B. D. Ripley, a quasinewton minimization method (BFGS) is used for the network's learning. We have taken the same nnet parameters as in [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF].

B. Evaluation of the neural network's outputs

A well-known problem, when using neural networks (or other regression methods), is overfitting: with enough parameters and enough training cycles, it is always possible to find a good fit for a given data set. So one may find a perfect fit for a chosen data sample, and then feel disappointed when the trained network makes wrong predictions on fresh data. So, we will systematically proceed as follows: train the network on a randomly chosen data sample (called train), then check the results, first on the same data sample, and second on a fresh data sample (called test) that was not used for the training.

In order to evaluate the outputs of several different models, we have to compare the neural networks predictions to the actual target values. We may use the fit criterion (cross entropy) but it does not reflect the influence of the number of weights (and biases) in the neural network. It is known (see [START_REF] Jordan | Neural Networks[END_REF]) that a network with too few weights may not be able to capture all the variations of the response to the input x, whereas a network with too many weights will more likely be subject to overfitting.In the next sections, we will compare several sets of input variables, of various sizes. Consequently, the number of weights in the network will not remain constant, and this variation will bias the results.

In [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF], we used the Akaike information criterion ( [START_REF] Greene | Econometrics analysis[END_REF]) to select the best model: AIC =2λ-2ln(L), where λ is the number of unadjusted parameters of the model (i.e. the number of weights and biases of the network), and ln(L) is the log-likelihood. However, fully connected neural networks are very greedy, as concerns the number of parameters: the number of weights and biases highly increases when adding new input variables. So we will rather use the Schwartz Bayesian information criterion (BIC = 2λ.ln(N )-2ln(L)), which assigns a higher penalty to the number of parameters in the model.

In our case, the BIC is written as follows, where N is the size of the data sample:

BIC =2λ.ln(N ) -2 N n=1 C k=1 t (n) k .ln(y (n) k ) (5)
In our previous work, we were able to compare predictions made on the train test and on the test set, which are of different size, by dividing the AIC by N , the number of data items: AIC avg = AIC N .In this paper, we will use BIC avg = BIC N . However, we must be aware that the correcting factor ln(N ) is different for both data sets, so it is meaningless to compare train and test data with this criterion. This is not of much importance, as we have already observed that the results on both data sets were highly consistent. In fact, we will make our model selection with the BIC criterion computed with the network's outputs predicted on test data only.

In addition to the numerical results provided by the Bayesian information criterion, we shall also consider the global proportion of correctly classified input vectors, and also the percentage of correct classifications for each class. One must be aware, however, that the rate of correct classifications is not the criterion being maximized by the neural network, so we should remain cautious when comparing the different classification rates. However, these percentages are easily understandable and may allow us to make some interesting statements on the results. Figure 2 shows, for train and test data, the evolution of the averaged Schartz Bayesian Information Criterion for several combinations of components. The mean values, averaged on five runs, are presented. This BIC curve confirms the similar results found in [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF] with the AIC.

V. Results

A. Components selection

The best trade-off between the model complexity (number of parameters in the neural network) and the benefit provided by the additional components is obtained with the {V ; C 1 ; ...; C 4 } model, that is: the sector volume, the "overal size" factor (aircraft count), the ground speed variance and vertical evolutions, the incoming flows, and the flow convergence and conflicts anticipation. Component C 5 (flow divergence) is either irrelevant, or redundant with C 4 (flow convergence and conflicts anticipation). Component C 6 (monitoring of aircraft vertical separation) is most probably irrelevant: in normal traffic situations, we may expect that all aircraft in close vertical proximity have previously been laterally separated by the controller. The same remark applies for horizontal proximity. This was already stated in [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]. So the aircraft in close proximity are not a significant factor in the decision to split or merge a sector. The anticipation of future conflicts may be more relevant.

This selection of the best model is very useful for a qualitative evaluation of the metrics. However, for practical purposes, it is not very convenient to compute the 27 complexity indicators, and then apply a transformation matrix to get the components. It would be useful to predict the sector status directly from a smaller subset of relevant complexity indicators.

B. Metrics selection

Let us now select the relevant metrics from the components of the best model {V ; C 1 ; ...; C 4 }. The initial variables subsets -the variables most correlated to each component axis -are detailed in table III of the PCA. We must be aware that the relative weight of a variable within a component (the linear correlation with the component axis) is not necessarily related to its relative influence on the quality of the neural network's prediction, as the relationship between the sector status and the explanatory variables may be non-linear. So we will again select the variables according to the Bayesian information criterion.

An iterative process is used to select the relevant variables. Each component is considered in turn. We have chosen Nb, the number of aircraft, to represent component C1, the overal size of the problem.

Starting with a basic set of variables {V, Nb},we build the best sequence of C 2 variables, by choosing at each step the variable which minimizes the BIC criterion. To do so, we consider each C 2 variable in turn (together with V and Nb), select the best one, add it to the sequence, and then consider the remaining variables, and so on until all variables have been tested.

The BIC curve is then drawn for this sequence, and the subset of variables corresponding to the minimum of this curve is added to the set of relevant variables. The process is then repeated with the next component, starting with the new set of relevant variables, and so on until all components have been considered.

The neural networks are trained on the train data, and the BIC is computed on test data. Five runs are made for each subset of variables. We will take the minimum value among these five runs, and use it in the selection process.

In the previous experiments ( [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF]) on components selection, the mean values and the minimum values lead to the same results. In the metrics selection, we have stated small differences between the two, probably due to the fact that some metrics belonging to a same component give fairly close results, like avg_vs, σgs gs , and σ 2 gs in component C 2 , for example. As the neural network with the best fit is closer to the "true" model than when considering the averaged value, we will now use the minimum value. C2)

The curve's minimum is reached with the combination {V, Nb, avg_vs}. This does not mean that the other C 2 variables are irrelevant: the results with {V, Nb, σ 2 gs }, and {V, Nb, σgs gs } were fairly close to the best result. So the ground speed variance and the ratio of standard deviation of speed to average speed are simply redundant with the average vertical speed avg_vs. This is also true for the numbers of climbing and descending aircraft, although they proved slightly less performant.

These results, showing the relevance of ground speed variance and altitude changes, are in contradiction with previous results presented in [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management[END_REF]. In this study, Masalonis, Callaham, and Wanke applied a logistic regression to subjective ratings and dynamic density metrics. The contribution of each metric was then assessed, for each ARTCC, by dropping the metric from the model. In particular, they found that the speed and altitude changes could be dropped without significantly affecting the model's predictions.

Our results contradict these statements, as we find that the component representing the speed variance and the altitude changes actually improves the BIC criterion. We have no ready explanation for this. It may be that, in [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management[END_REF], the sector volume alone captures the fact that some sectors are dedicated to a specific type of traffic: pre-approach sectors with evolutive aircraft are generally smaller than pure en-route sectors. As we use sectors or groups of sectors of various sizes and statuses (merged, armed, or split), there is no correlation between the sector volume and the nature of the traffic within the sector. So the speed variance and the altitude changes become relevant in that context. But, as we don't know the particulars of [START_REF] Masalonis | Dynamic density and complexity metrics for realtime traffic flow management[END_REF], it is difficult to be sure.

Let us now consider the variables of component C 3 , applying the same procedure, and starting from the new set of relevant variables {V, Nb, avg_vs}. Figure 4 shows the evolution of the BIC avg for the best sequence found for component C 3 (see table V). The minimum is reached for the subset {V, Nb, avg_vs,F 60 ,F 15 }.

Subset Indicators C3_a

V , Nb, avg_vs, F60 C3_b

V , Nb, avg_vs, F60, F15 C3_c

V , Nb, avg_vs, F60, F15, F5 C3_d

V , Nb, avg_vs, F60, F15, F5, F30 C3)

The incoming flow with a time horizon of 15 minutes belongs to the subset of relevant variables. This is consistent with the fact that strips are sent to the controller's working position about 15 minutes before the aircraft actually enters the sector.

We were rather surprised to see that the incoming flow of traffic entering the sector in the next hour also belongs to the set of relevant variables. This incoming flow is probably less linked to the variations of the number of aircraft, and it may somehow smooth these variations.

The incoming flow within the next 5 minutes does not provide any improvement. The first step of the BIC minimization process showed that the BIC avg criterion increased when adding F 5 to the subset {V, Nb, avg_vs}, whereas all other flows improved the criterion. So it is either highly redundant with the previous variables, or irrelevant. An interesting result is that the metrics related to the anticipation of conflicts, for converging aircraft in the horizontal plane, do not bring any supplementary benefit, when added to inter_hori, the number of potential crossings with angle greater than 20 degrees. Notice that all these metrics -anticipation of conflicts and flow convergence -were already put together in a same component by the PCA. This leads to think that the main feature captured by this set of metrics is the fact that the pairs of aircraft within the sector are in a converging configuration, or not.

We have a feeling that other conflict detection indicators (not only for converging flights, and not only in the horizontal plane) may be useful, and may be more orthogonal to the indicators we have implemented so far. This may be the subject of future work.

To conclude on the metrics selection, we end up with a subset of only 6 most relevant indicators {V, Nb, avg_vs,F 60 ,F 15 , inter_hori} among the 27 metrics, plus the sector volume, that were considered in this study. This allows to establish a direct relationship (equation IV-A, with the weights and biases of the trained network) between these six indicators and the probability to belong to a sector status class.

C. Classification rates

VI. Conclusion

In conclusion, the proposed method, using neural networks and sector status records, allows selecting the most relevant air traffic complexity metrics, among several indicators proposed in the literature.

In a previous work ( [START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF]), a principal component analysis reduced the scope of the study to 6 main components , starting from 27 metrics to which we added the sector volume. It was then found, by applying neural networks to the principal components and to the sector statuses, that only the first four were significantly related to the decisions to split or merge sectors.

In the current paper, we focused on the subset of metrics belonging to the four relevant components. Still using neural networks and sector statuses, we proposed a selection process, minimizing the Bayesian information criterion, which allowed to select a small subset of 6 relevant metrics: the sector volume V , the number of aircraft within the sector Nb, the average vertical speed avg_vs, the incoming flows with time horizons of 15 minutes and 60 minutes (F 15 , F 60 ), and the number of potential crossings with an angle greater than 20 degrees (inter_hori). Let us note that the other metrics may be either irrelevant, or simply redundant with the most relevant indicators.

We established a fairly simple equation (equation IV-A with the weights of the trained network), allowing to predict the sector status (merged, armed, or split) from these relevant metrics. The rates of correct predictions were above 83% (global rate), and highly consistent when considering train and test data.

The neural network approach used in this study seems appropriate for the granularity we are interested in and the foreseen applications, either strategical (sector design) or pre-tactical (sector planning). We are also fairly confident that decisions to split or merge sectors may allow to assess the instantaneous workload as well, and could therefore be used to improve tactical tools (PRESAGE). To this end, other statistical methods should be investigated to take into account the serial correlation of sector status, looking closely at the sector splitting times. We plan to tackle this issue in a close future, using dynamic discrete choice models.

Finally, another issue that we intend to address, in relation to the complexity indicators, is the prediction of optimal sector configurations. Previous works ( [START_REF] Gianazza | Optimization of air traffic control sector configurations using tree search methods and genetic algorithms[END_REF], [START_REF] Gianazza | Optimal combinations of air traffic control sectors using classical and stochastic methods[END_REF]) proposed several algorithms to compute optimal sector configurations, using sector capacities and incoming flows. The output of the neural network is a triple of probabilities, allowing to decide when a sector should be split, or merged. We may derive a realistic workload indicator -and also threshold values -from these probabilities, which could be used to compute optimal sector configurations.

Metrics inspired from the CREED project

The work of P. Averty on conflict detection [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] inspired a set of indicators. One of the ideas in [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] is that conflict perception is "plannar". The author thus defines for converging pairs of aircraft the following quantities

• Ed : minimum horizontal distance between aircraft.

• Efl : horizontal distance when the aircraft are vertically separated (after the crossing). • Da : the "anticipation degree", i.e. the distance between the faster aircraft and the intersection of the aircraft trajectories (in the horizontal plan). We replace this variable to a modified Da, DaC, which stands for the greater distance between one of the aircraft and the point where, horizontally, the distance between aircraft is the smallest. For explanations about this substitution, we refer to [START_REF] Guittet | Analyse descriptive des indicateurs de complexité du trafic aérien à partir des données image et courage[END_REF].

Originally, these quantities are defined to describe conflict perception. To translate the idea of [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] in terms of traffic complexity, we assume that a conflict is all the more critic that the expected separation (Ed and Efl) and the anticipation (DaC) are small. We thus set

creed = 1 α Da +(1-α)(β Ed +(1-β) Efl) , (18)
where α and β are parameters in [0; 1] 4 . Finally, aircraft pairs considered in [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] are such that vertical separation occurs prior to separation, as the converse situation is avoided as much as possible by controllers. Accordingly, the complexity associated with these latter pairs is likely to be greater and we distinguished the two kind of conflicts by summing the quantity fintroduced in [START_REF] Kopardekar | Dynamic density: A review of proposed variables[END_REF] on both sets of aircraft, thus creating two distinct indicators, creed_ok ("good pairs") and creed_pb ("bad pairs").

25th Digital Avionics Systems Conference October 15, 2006 4 As for now, these parameters are set equal to 0.5, but are meant to be adjusted and possibly vary with DaC to reflect the results of ongoing research [START_REF] Averty | Work in progress, presented at an internal SDER seminar[END_REF].

Figure 1 .

 1 Figure 1. Example of a feed-forward network with one hidden layer

Figure 2 .

 2 Figure 2. Values of the BICavg criterion for the different sets of components

Figure 3 .Figure 3

 33 Figure 3. Evolution of the BICavg criterion for component C2, for the sequence avg_vs, σgs gs ,σ 2 gs ,N cl ,N des

Figure 4 .

 4 Figure 4. Evolution of the BICavg criterion for component C3

TABLE II NUMBER

 II OF MEASURES AVAILABLE, ON THE 1ST JUNE OF 2003.

  2 gs ,N cl ,N des . The x axis labels on this figure are detailed in table IV.

	Subset Indicators
	C2_a	V , Nb, avg_vs
	C2_b C2_c C2_d	V , Nb, avg_vs, V , Nb, avg_vs, V , Nb, avg_vs,	σgs gs σgs gs , σ 2 gs σgs gs , σ 2 gs , N cl
	C2_e	V , Nb, avg_vs,

σgs gs , σ 2 gs , N cl , N ds

TABLE IV SUBSETS

 IV OF THE SEQUENCE SHOWN ON FIGURE 3(COMP.

TABLE V SUBSETS

 V OF THE SEQUENCE SHOWN ON FIGURE 4(COMP.

  Figure 5. Evolution of the BICavg criterion for component C4 Figure 5 and tableVIshow the results for component C 4 . The minimum for this component is reached when adding inter_hori. The additionnal variables of the rest of the best sequence increase the BIC criterion. However, the other variables also performed quite well, when put in first position in the sequence. This means that they are not irrelevant, but redundant with inter_hori.

	Subset Indicators
	C4a	V , Nb, avg_vs, F60, F15, inter_hori
	C4 b	V , Nb, avg_vs, F60, F15, inter_hori, creed_pb
	C4c	V , Nb, avg_vs, F60, F15, inter_hori, creed_pb,
		creed_ok
	C4 d	V , Nb, avg_vs, F60, F15, inter_hori, creed_pb,
		creed_ok, Conv
	C4e	V , Nb, avg_vs, F60, F15, inter_hori, creed_pb,
		creed_ok, Conv, insen_c
			0.85
			0.84
		BIC avg	0.83
			0.82
			0.81
		V+Nb+avg_vs+F60+F15	C4_b	C4_c	C4_d	C4_e
				Metrics subset

  TableVIIshows the correct classification rates obtained, on train and test data, with the subset of most relevant indicators. The results with test data are rather consistent with the ones obtained with train data, as was already the case when studying the components ([START_REF] Gianazza | Evaluation of air traffic complexity metrics using neural networks ans sector status[END_REF]).

	Set	Global	Merged	Armed	Split
	Train	83.47%	89.07% 65.58% 91.91%
	Test	83.24%	89.14% 64.56% 91.75%

TABLE VII CORRECT

 VII CLASSIFICATION RATES FOR THE {V, Nb, avg_vs,F60,F15, inter_hori} MODEL, WITH A

	I6H15O3 NETWORK

Reasonable weights were given by P. Averty and M. Tognoni.

APPENDIX A: Complexity metrics

Delahaye and Puechmorel metrics

To present the geometrical indicators introduced in [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], we need to define several quantities:

• The vector representing the distance between two aircraft is denoted by ---→ X i X j where X i (resp. X j ) stands for the location of aircraft i (resp. j).

• The "oblical" distance between two aircraft (i and j) is denoted by

where <. , .>stands for the appropriate scalar product.

between two aircraft. • The derivative of the "oblical" distance between two aircraft is denoted by v ij and writes

• We introduce a weighting function f . As suggested in [START_REF] Chatton | Etudes de nouvelles métriques de complexité de la circulation aérienne[END_REF], we used

with α =0.002, β =0.01, the d ob being expressed in nautical miles. These indicators are defined pointwise. To get a value on the controlled airspace, they have to be averaged on the different aircraft. In [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], a density indicator is defined as follows

Two indicators are introduced to reflect the variability in headings (track_disorder) and speed (speed_disorder).

There are defined as

Indicators Div et Conv respectively describe convergency and divergency of the aircraft in the controlled sector.

Indicators Sd + and Sd -are designed to set a weight on potential conflicts that are difficult to solve. These "sensitivity" indicators are defined by

Note that components of the gradient are weighted so as to reflect the difficulty of the respective manoeuvers 3 .As observed in [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], a situation with a high "sensitivity" is easier to resolve for the air controller than one with a low "sensitivity". As these indicators "increase" with the number of aircraft, it is unclear whether they actually are "complexity" or "simplicity" indicators. We thus define a last pair of indicators, insen_c and insen_d as

Modified PRU metrics

The work conducted by SDER-RFM for the Performance Research Unit (citeRFM), though initially designed to compare ATC centers on a daily basis, inspired the following indicators :

• inter_hori: number of potential crossings (irrespective of the aircraft direction on their trajectories) with angle greater than 20 degrees.

• avg_vs: this is simply the average vertical speed of controlled aircraft.


S Puechmorel

4D TRAJECTORIES: A FUNCTIONAL DATA PERSPECTIVE

3D or 4D trajectories are fundamental objects within the frame of ATM and can be defined as mappings from a bounded interval to the space 3 R . However, this functional aspect is almost never used in applications. A great improvement over existing procedures for statistical analysis of trajectories can be obtained by explicitly considering trajectory data as functional data. Since raw functional data is by essence infinite dimensional, no computation can be made unless a finite representation has been found: all the complexity of functional data analysis is hidden in this stage. The present paper will first analyze the functional aspect of trajectories in order to extract operationally relevant quantities then will give a way of representing optimally aircraft trajectories with a finite number of parameters, allowing further processing. Applications to trajectory prediction and classification will be briefly discussed at the end.

Introduction

Aircraft trajectory is one of the most fundamental objects within the frame of ATM. However, partly due to the fact that aircraft positions are most of the time represented as radar plots, the time dependence is generally overlooked so that many trajectory statistics conducted in ATM are spatial only. Even in the most favorable setting, with time explicitly taken into account, trajectory data is expressed as an ordered list of plots labeled with a time stamp, forgetting the underlying aircraft dynamics. On the opposite, functional data analysis will work with the whole trajectory as a mapping from a time interval to the state space (most of the time 3 R or 6 R when speed is considered as being part of the aircraft state). There is quite a long history of functional data in statistics [1,[START_REF] Ramsay | Functional Data Analysis[END_REF], but to our knowledge no application to the ATM field has been done yet. Furthermore, this way of thinking induces a change of paradigm encompassing more than bare statistics: functional representation makes an assumption on the underlying phenomenon giving rise to the observed data so that a great wealth of information can be taken into account. As an example, aircraft trajectories are mainly made of straight lines and arcs of circles. Using a linearized model of flight dynamics, such trajectories can be encoded using only piecewise constant mappings (namely the commands applied to the system) so that the whole data can be expressed using only a small number of samples. On the opposite, the collection of radar plots describing the same trajectory can have tenths more samples, nearly all of them redundant. As a consequence, data analysis based on the first representation will use only the salient features of the trajectory while with the second a lot of useless information will be injected.

In the present paper, we will analyze first the possible representation for functional data, then some application will be presented.

Air Traffic Management Applications

Regular aircraft trajectory representation (ordered list of plots labeled with time) does not contain typical ATM operational information due to the limited information given by a list of plots. The new representation described in this paper takes into account typical ATM operational features such as curvature and torsion. These relevant features are naturally encoded in this new representation and can be used for turns detection or any other maneuvers. Based on these new trajectories representation many ATM operational applications may be investigated.

Aircraft Trajectories Data Compression

As it has been previously mentioned, ATM systems manage aircraft trajectories and control them in order guaranty safety and airspace capacity. Currently those trajectories are represented by the mean of plot lists which are manipulated by ATM software. Every day, all aircraft trajectories are registered into large database for which huge capacity is needed. Based on this new trajectory 1-4244-1108-4/07/$25.00 ©2007 IEEE.

1.C.6-1 representation for which redundancy has been removed, the trajectories database may be strongly improved from the capacity point of view. As a matter of fact, only spline decomposition coefficients have to be registered in the database then reducing the needed capacity. This compressed trajectory format may also be used for improving the trajectories transmission between ATM entities.

Aircraft Trajectories Distance Computation

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. This new trajectory representation will enable the definition and fast computation of a rigorous mathematical distance between aircraft trajectories. Based on this new distance computation several ATM applications may be implemented.

Aircraft Model Inference

All aircraft models are based on ODEs(Ordinary Differental Equation), including tabular ones. Control input includes condition and model parameters. The model refinement (and computational complexity) ranges from tabular to many degrees of freedom. 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. In order to find the closest model in this trajectory space, a reliable trajectory distance is needed. The model inference problem has to solve the accuracy-smoothness 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. All deformations from the goal trajectory to the trajectory set are taken into account but many of those transformations are not feasible from a physical point of view. The solution consists in defining a smooth path between trajectories by an infinitesimal displacement field and by adding a field smoothness requirement to the energy functional. On the other hand, if we decide to control the maximum distance between the goal trajectory and a given parametrized model, it is possible to adjust the most economical model to the current phase of flight. As a matter of fact, computation effort is proportional to the degrees of freedom used and to the evolution equation used. Then, the relevant question is the following: Can we switch from different level of accuracy depending on flight phases and computational power available? The solution consists in switching from one model of given accuracy to another by following a shortest path in the trajectory space.

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 | Us Request Evaluation tool (URET) Conflict Pro 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 T MA[END_REF], controller advisories for metering and sequencing [START_REF] Swensen | Design and Operational Evaluation of the Traffic Management Advisor at the Forth Control Center[END_REF][START_REF] Coppenbarger | Design and Development of the EnRoute Descent Advisor (EDA) f Arrival Metering[END_REF],traffic load forecasting [START_REF] Meckiff | The Tactical Load Smoother for Multi-Sector Planning[END_REF][START_REF] Masalonis | Using Probabilis Demand Prediction for Traffic Flow M 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 Coord Tools[END_REF][START_REF] Evans | Reducing Severe Weathe Delays in Congested Airspace with Weather Support for Tactical Air Traffic Budapest Hungary[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 continue as active areas of research and development (e.g. [14][15][START_REF] Mondoloni | Trajectory Modeling Accuracy for Ai Traffic Management Decision Sup Toronto ON[END_REF][START_REF] Ryan | Review of Trajectory Accuracy Methodo 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. Since we may prescribe smoothness accuracy of a simplified model relative to a finer one, may be computed. Finally such trajectory distance allows comparison between TP models. 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

1.C.6-2
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 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 (biallocation). 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 over a long period of time in a dense area, it is very easy to identify 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 split 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).

Trajectories as Mappings

In the following all aircraft trajectories will be described as mappings from a time interval [] ab to a state space with either

EE 3 R or 6 R
depending on the fact that speed is assumed to be part of aircraft state or not. Extension to trajectories on a sphere (typically long haul flights) will be sketched only.

Notations and Terminology

The reference for this section is [START_REF] Berger | Differential geometry: Manifolds, Curves and Surfaces[END_REF]. Let [] ab E be a trajectory. The origin of the trajectory is () a and the destination is () b .

Those two points are called the endpoints of the trajectory. All trajectories are assumed to be at least continuously differentiable (class ) so that the length of a trajectory

R since it is not 1.C.6-3
possible for an aircraft to stop except at the endpoints of the trajectory.

Remark 3. While the case

is very natural and intuitive, care must be taken when since all the preceding definitions apply in a completely different setting: for example, the non singularity assumption does not implies nowhere zero speed, but only that speed and acceleration cannot both vanish at the same time.

() () Ks s

The curvature can be explicitly computed even if the curve is not parametrized by arclength. The general formula is:

3 () () () () tt Kt t
with the vector cross product. Curvature is of primary importance for ATM related studies since as mentioned before aircraft trajectories are mainly made of straight lines and arcs of circle and so have piecewise constant curvature. If at point t the curvature is not zero, the curve is said to be biregular at t. For a curve parametrized by arclength, the unit normal vector s is defined at all biregular points by: () () ()

s

sT s s

The real number is called the torsion of the curve at s and represents an obstruction for the curve to be planar. As for the curvature, it is possible to compute the torsion even if the curve is not parametrized by arclength:

Ts 2 () ( () tt tt (( ) ) ) () ()

det t Tt

Torsion is not so useful as curvature for enroute data analysis since only a few number of trajectories have non zero torsion. However, it is very relevant in terminal areas.

Remark 5. The case is again very different, since the geometric meaning of curvature and torsion is not obvious in this setting. Furthermore, the extra degrees of freedom will impose using higher order derivatives in order to build up an equivalent description. A complete treatment goes beyond the scope of the present paper and has little interest for our purpose (in practical applications, the speed information, when available, is used to improve estimates of curvature and torsion and not to study a trajectory in

6 ER 6 R ).
The fact that the last two quantities are nearly piecewise constant functions for aircraft trajectories makes the previous representation very efficient for ATM applications. Most of the time, the relevant assumption will be that the torsion is piecewise the first three derivatives of the trajectory are available. In applications however this is nearly never the case since even the speed has to be computed. Moreover, real data are corrupted b noise, making evaluation of derivatives quite cumbersome. To obtain good estimates of curv and torsion, we will first compute an expansion of the trajectory on an infinite dimensional basis that will be truncated to obtain a finite representation of the trajectory. Since we know that the relevant information is hidden behind curvature and tors that is quantities related to first, second and third derivative, the decomposition basis has to be chos so that those values are easy to compute. t ( 0-th order derivative is the polynomi elf d so this conditions can be used to explicitly take the value N the al its an y into account). Although it is possible to compute an interpolating spline based on the previous derivations, it is more enlightening to take a different approach for trajectories representation. First, one can note that a 0-th order spline is indeed a piecewise constant function and thus is entirely defined by the interpolation data (no smoothness is required). It is then clear that the interpolating spline can be written as a weighted combination of elementary piecewise constant functions:

1 1 , 0 1 ii N i tt i

St y t

In the general case, we will seek for a deco inimally mposition on basis functions that are m compactly supported splines. For regularity d and a spline function made from N polynomials, smoothness requirement will induce

R .

The optimal has expansion: nstant curvature and torsion, the choice of L will be made in order to satisfy as well possible this requirement. Ideally, L should b derivative of the curvature (or the torsion), but is nonlinear, an optimal L cannot be found. The standard way of dealing with the problem is to require only that the square norm of the third derivative be as small as possible, thus giving: e components: all the results previously given will apply. Note that many standard spline smooth algorithms use instead cubic splines that corres to an approximate minimum total curvat Because of the approximation made on the criterion, curves reconstructed by quintic splines will not be of piecewise constant curvature. In order to have a better approximation, it is convenien add a penalty term so that the norm of the ve be as constant as possible: in this case, curves tend to be parametrized by arclength, thus giving a better approximation of curvature by the second derivative. The corresponding composite criterion is: 00 00 D LD 00D 0 al with intensively studied in t case)

a positive real number controlling the importance paid to consta velocity. Functions optim respect to this new criterion have been he literature (at least for the nt cubic and are named splines under tension ( is the tension parameter). Increasing results in curve with sharper turns and more straight line segments, while setting 0 a gives the initial quintic (or cubic) splines. As it has been pointed out by K. Mehadehbi, it is interesting to make depend on time, so that it is possible to more accurately model turns by relaxing tension while low curvature areas will be associated to high tension, yield almost straight lines which is operationally sound. The choice of the tension parameter remains anyway a highly heuristic process. In order to more accurately take into account the real curvature and more generall y the 1.C.6-8 shape of the trajectory to be reconstructed, a different approach combining a shape optimiz an a velocity interpolation will be taken. Recal shape parameters like curvature or torsion are easil defined only when the curve is parametrized b arclength. Assuming that such a parametrization has been obtained, the interpolation criterion has to be modified: with time parametrization tha curves that ca be obtained by rotations of th same curve will be equally penalized), it can be g clu Schur lemma that PP has to be ing the use of vector valued kernels. Remark 6. Since most of the relevant P for aircraft trajectory studies are diagonal, this restriction has little impact: interpolation can be done with scalar kernels component-wise. t d ermore, computational speed is greatly increased when using such kind of kernels: a up of 9 is theoretically possible since the complexity of linear algebra algorithms use the coefficients of spline expansion is 3 oN with N the dimension of the matrices involved. With general vector kernels, this amounts to a complexity of order 3 32 7 while the diagonal cas es to 3 independent resolutions of one dimensional problems, thus yielding an overall complexity of order 3.

In [START_REF] Smola | Regularized Principal Manifolds[END_REF], it is proposed to find a solution to the principal curve problem by assuming that all curv e reduc es have an expansion as a weighted sum of a reproducing kernel: mapping: it will be obtained again using a spline interpolation, but this time using a well known cubic spline interpolation. It must be noted that monotonic splines exist, thus insuring that the parametrization has by construction the required properties, however, on practical cases, it appears that the simple cubic splines algorithm produces admissible results with a much lower computational cost. For that reason, only this last option has be considered in our final implementation of spline smoothing.

Using Functional Data

Once a tractable (i.e. finite dimensional) representation of the original function classical vector algorithms to th framework. For the purpose of aircraft trajectory analysis, two applications have been investigated

Trajectory smoothing and prediction. Trajectory classification.

We will now briefly explain how the use functional data representation can improve o existing procedures. Trajectory s n of sampled data is obtained by solving a smoothing problem. For prediction, sort term forecasting can be again directly obtained by evaluating the functional representation at t desired time. However, smoothing splines, eit standard or shape based are generally poor predictors for medium to long term forecasting: some knowledge about the future behavior of the trajectory has to be incorporated to the model. Some information from flight plan can be added is available by considering those intended pos measures, but with a low confidence level (this is obtained simply by weighting the correspond error term in the smoothing problem). Since 4D posi are generally badly defined in flight plans, it is far better to use the regularized principal curve algorithm first, then to use the 1D cubic spline smoothing to infer the time/arclength mapping, with a low confidence on values obtained from the flight plan). For long term forecasting, only statistical procedures are relevant. Since it amounts to f the most probable trajectory based on observatio it is in fact based on the next item that is trajectories classification.

Classification of trajectories is a statistical procedure aimed at splitting a dataset into subsets with high internal coherence and at finding a best representing member for each of those subsets. Many algorithm n one of those two classes:

Supervised classification, where the number of classes is known and an already classified learning dataset is available.

Unsupervised classification where the classes are inferred by the algorithm.

Supervised classific applied in the functional context since almost algorithms require only a measure of similari In our setting, such a measure can be tional data (for operational purpose, an interesting distance is obtained by summing the 2 L norm of the differences of trajectories, deriva of trajectories and second derivatives of trajectories). For the unsupervised case, most algorithms include a heuristic step: we have investigated such algorithms for the purpose of craft trajectory analysis.

Conclusion and Future Work

Because trajectories are functional objects i highly desirable to stay with logy in the frame of ATM are nearly all the time based on sampled data, a way of optimally inferring a functional model has been presented. Since relevant operational quantities are related to the shape of trajectories, it is interesting to first model the data as a regularized principal curve, then to fit an ad-hoc parametrization to match the 4D aspec This algorithm improves over the smoothing spline or splines in tension algorithms already in use. Furthermore, it is much easier to enforce piecewise constant curvature or torsion using shape models. The price that has to be paid for this extra flexibility is a more computationally intensive procedure, involving an iterative optimization stage. Taking into account the computation power of nowadays computers, it turns out that computational time penalty compared with standard smoothing spline is acceptable. A extension of the shape algorithm using an optimal tension parameter is currently under development. At the same time, a dataset of arrival trajectories at Paris CDG will be analyzed using the functional data approach.

3 ER 6 ER 3 ER

 363 The arclength parametrization allows to define very important geometrical quantities when .

  1.C.6-4 constant, with value 0 in all flight phases except those occurring in terminal areas. Curves with constant torsion and curvature are helices, so th our characterization of aircraft trajectories by piecewise constant curvature and torsion curve induce a representation by pieces of helices.All the previous derivations rely on the fa at

1 N

 1 A nice introduction thing applied to aircraft trajectories can b found in a technical report from K. Mehadhebi [2 In their simplest form, splines are piecewise polynomial approximations to real function, w the extra regularity assumption that the reconstructed mapping has to be of presc class r C . In the case of functions taking their va in p R , splines will be constructed coordinate-wise. Let , I ab be a real interval and let: points are called the knots of the spline). A polynomial spline on the subdivision is a mapping

1 R

 1 obtain a RKHS for our problem, we will assume that all trajectories are defined on the s whose derivatives up to order m exist (in the weak sense) and are square summable. Let L be a differential operator of order m. The n: the previous sense). Furthermore, the Green function of L gives immediately access to the reproducing kerne with R l of . Recalling that our purpose is to obtain e co as e th unfortunately, since the expression of the curvature with D the operator of differentiation with respect to time. Simple computation shows that this particular choice gives piecewise polynomial splines of degree 5 on each of the three

  One must be careful with the respective definitions of trajectories and curves: a curve is defined on an open interval and thus has no endpoints. Nevertheless, any trajectory

						C	1
						[] ab	E	is well defined
	as:
				() lt dt () b a
	If	() 0 t	for some	( b ) ta	the point t is
	said to be singular. A parametrized curve of class
	(or more concisely a CC p	p	curve) will be a C	p
	mapping from an open time interval	() ab to the
	state space	E	with no singular points. Any	C	1
	curve can be parametrized by arclength. Let
	() ab		E	be such a curve. Defining the
	mapping ()( 0 bl () ) sa	by:
				() st	t a	() t dt
	we see that by the non singularity assumption on ,
	() st	() 0 t	for any ta, so that is ( ) b s
	an invertible mapping. Now,	s	1	is a mapping
	from the open interval	(0 l	() )	to E satisfying:
		11 () ( ss )(	s	1 )1
	In the following, we will simply write
	() ss	( 0 ( ) ) l
	:, ab		E has an associated curve, namely
	:, ab		E . It is generally more convenient to
	deal with curves to avoid special treatment of the
	endpoints.
	Remark 2. The non singularity assumption on
	the underlying curve is very natural when dealing
	with aircraft trajectories in 3

for a curve parametrized by arclength, dropping the variable t . Remark 1.