Roman Goulard

Anna Verbe

Jean-Louis Vercher

Stéphane Viollet
email: stephane.viollet@univ-amu.fr

Role of the light source position in freely falling hoverflies´stabilization performances

Keywords: flight stabilization, insect vision, Dorsal Light Response, hoverflies, dipterans

The stabilization of plummeting hoverflies was filmed and analysed in terms of their wingbeat initiation times as well as the crash and stabilization rates. The flies experienced nearweightlessness for a period of time which depended on their ability to counteract the free fall by triggering their wingbeats. In this paradigm, hoverflies' flight stabilization strategies were investigated here for the first time under two different positions of the light source (overhead or bottom lighting). The crash rates were higher in bottom lighting conditions than with top lighting. In addition, adding a texture to the walls reduced the crash rates only in the overhead lighting condition. The position of the lighting also significantly affected both the stabilization rates and the time taken by the flies to stabilize, which decreased and increased under bottom lighting conditions, respectively, whereas textured walls increased the stabilization rates under both lighting conditions. These results support the idea that flies may mainly base their flight control strategy on visual cues and particularly, that the light distribution in the visual field may provide reliable, efficient cues for estimating their orientation with respect to an allocentric reference frame. In addition, the finding that the hoverflies' optic flow-based motion detection ability is affected by the position of the light source in their visual field suggests the occurrence of interactions between movement perception and this visual vertical perception process.

Introduction

Several flying insects including dipterans are known to keep their dorsal surface oriented upwards by holding the brightest part of their environment, which usually shines from above, in a constant position in their visual field [START_REF] Mittelstaedt | Physiologie des Gleichgewichtssinnes bei Fliegenden Libellen[END_REF]. This reflex, which is known as the Dorsal Light Response (DLR) has also been described in detail in fish [START_REF] Fraenkel | The orientation of animals: Kineses, taxes and compass reactions[END_REF]. The importance of the orientation of an artificial horizon in blowflies' head roll orientation processes has also been previously established, as well as the fact that these insects probably do not use gravity information to perform this task [START_REF] Hengstenberg | Localization and Orientation in Biology and Engineering[END_REF], which has been assessed in freely flying hoverflies based on free fall experiments [START_REF] Goulard | To crash or not to crash: how do hoverflies cope with free-fall situations and weightlessness[END_REF], contrary to the well-known negative gravitaxy behaviour observed in walking Drosophila [START_REF] Kamikouchi | The neural basis of Drosophila gravity-sensing and hearing[END_REF][START_REF] Hidehiko K Inagaki | Methods for quantifying simple gravity sensing in drosophila melanogaster[END_REF]. These results suggest that visual processes predominate over gravity-based ones in the strategies used by flying flies to stabilize their flight, and support the idea that there exists some kind of vertical reference frame in flies' brains based on the DLR. However, although this reflex had been found to be closely linked to the head roll steering mechanism, this situation has been established only in tethered Calliphora and Episyrphus [START_REF] Hengstenberg | Visual Motion and its Role in the Stabilization of Gaze[END_REF][START_REF] Goulard | Behavioural evidence for a visual and proprioceptive control of head roll in hoverflies (episyrphus balteatus)[END_REF] and has not been studied so far during free flight.

In the present study, it was therefore proposed to assess the effects of the change in the light source position on freely flying dipterans' stabilization performances. These performances were tested using a free fall procedure under four different visual conditions in which two differently textured lateral walls (uniform and textured) were combined with two different lighting sources (overhead or bottom lighting). The position of the light source was found to be a crucial factor for hoverflies to be able to regain a suitable flight attitude after a free fall. In addition, the results obtained in this study support the idea that the use of lateral visual cues such as 2-D chequerboard patterns generating optic flow (OF) may be involved in hoverflies' attitude and lift control processes [START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF][START_REF] Goulard | Modeling visual-based pitch, lift and speed control strategies in hoverflies[END_REF]. This study clearly shows not only that the light gradient perception and OF-based control processes are interlinked, but also that the static cues consisting of the light source position predominate over the insects' OF-based control processes.

Methods

Animals

Hoverfly pupae (Episyrphus balteatus) were purchased from Katz Biotech AG, Baruth, Germany. To magnetically maintain the animals in resting position (see in figure 1A), a piece of entomological pin approximately 5 mm long was glued to the dorsal part of the animals32 thorax, perpendicularly to their longitudinal axis: the pin (≈ 5mg) weighed approximately 15% of the hoverfly's mass (≈ 35mg). The insects' flight ability was checked in the breeding cages throughout the experiments. 39 hoverflies (19 in Halogen/LED experiments and 20 in the control experiments) aged from 3 to 28 days were tested (9 males and 10 females in the Halogen/LED experiments and 2 males and 18 females in the control experiments).

Experimental procedure

Hoverflies were subjected to free fall conditions in a modified version of the setup previously presented in [START_REF] Goulard | To crash or not to crash: how do hoverflies cope with free-fall situations and weightlessness[END_REF]. In the present version, the box was illuminated alternately from above (Top Lighting: TL) with a white halogen light (Kaiser Studiolight H) and from below (Bottom Lighting: BL) with a white backlit LED panel (Phlox, 50x50cm 2 ) featuring a uniformity as high as 95% and two peaks (450nm and one at around 550nm) that match the spectral sensitivity of the hoverfly's (Erisalis tenax ) photoreceptor cells [START_REF] Ga Horridge | Fly photoreceptors-ii. spectral and polarized light sensitivity in the drone fly eristalis[END_REF]). in size). The irradiance was measured in both illumination conditions with an ILT1700 radiometer (International Light Technologies) under both experimental conditions (with textured walls: CT and CB conditions) by orienting the light probe (SED033, visual field 3°) of the radiometer towards either the illuminated side or the opposite side. The irradiance measured in the CT condition (direct/indirect measurements) was 1.12.10 -8 /3.15.10-9 W.cm -2 and 1.10.10 -9 /2.33.10 -10 W.cm -2 in the CB condition with the LED light. In the control experiments with two halogen lights, the irradiance was set at the same value (direct measurements) (1.24.10 -8 W.cm -2 ).

A total number of 262 falls were conducted among the four different conditions (see figure 1B), and 91 additional falls were conducted in the control experiments. At each experimental session, a hoverfly was exposed to the four environments consecutively in random order.

Each hoverfly could undergo several experimental sessions, but no more than once a day in order to prevent the occurrence of any habituation or fatigue effects. We always checked between experimental sessions whether the hoverflies equipped with their glued pin were able to fly in the breeding cages.

Image analysis

The horizontal and vertical 2-D positions of the hoverflies' centre of mass moving over a uniform background were recorded using a custom-made image-processing program running under MATLAB. The fly's speed was calculated from the positions recorded by applying a Savitzky-Golay procedure (order 2, window: 51). Stabilization was determined automatically when the fly reached a positive vertical speed without touching either a wall or the ground.

Statistical analysis

Data were analysed statistically using a generalized linear mixed-effects model procedure ('glmer' in R v3.2.3) and selected using the Akaike Information Criterion ['AIC' [START_REF] Akaike | A New Look at the Statistical Model Identification[END_REF]].

Results

As observed previously [START_REF] Goulard | To crash or not to crash: how do hoverflies cope with free-fall situations and weightlessness[END_REF], hoverflies subjected to free falls initiated their flight after approximately 100ms in both uniform and contrasting wall texture conditions (p = 0.12; F = 2.4727; Figure 2A). However, top lighting conditions significantly decreased the reaction wingbeat triggering) times (p < 0.01; F = 9.4352). The flies' performances in the CB condition differed significantly from those observed in the CT (p < 0.05; z = -2.863) and UT (p < 0.01; z = -3.294) conditions, but the differences in the mean times did not exceed 20ms ( ∆WB (ms): CB = 131.2500; CT = 107.7500; U B = 116.5152; U T = 104.3087). It is worth noting that during our control experiments, the wingbeat triggering times (see Figure 2A) were significantly shorter by around 20ms (p < 0.001; F = 19.6268), but that this did not significantly reduce the difference in the effects observed between CB and CT (p = 0.39; F = 0.7303).

Bottom lighting conditions induced a much larger number of touchdowns on the floor (Fig- ure 2B), amounting to approximately 60% of all the trials, than under overhead lighting conditions (p < 0.001; Chi 2 = 35.8369), which enabled the hoverflies to avoid crashing in 70-75% (UT) to 90% (CT) of the flights. It is worth noting that the crash rates were quite similar between Halogen-LED and control experiments (Effect of LED light: p = 0.68; Chi 2 = 0.1694; interaction with light position effect: p = 0.35; Chi 2 = 0.8571) which confirms the validity of using a LED panel to stimulate light-dependent stabilization behaviour. In addition, the crash rates were not significantly affected by the presence of textured walls (p = 0.18; Chi 2 = 1.7578), whereas a significant interaction was found to occur between the lighting and texture conditions (p = 0.01; Chi 2 = 6.5524). The presence of a 2-D chequerboard pattern on the walls significantly decreased the crash rates under overhead lighting conditions (Post-hoc Tukey contrast, CT vs UT: p < 0.05; z = 2.838) but not under bottom lighting conditions (Post-hoc Tukey contrast, CB vs UB: p = 0.95; z = -0.525).

In the subsequent analysis, stabilized flight was taken to occur whenever the fly adopted a positive vertical speed, corresponding to a rising flight, without subsequently crashing onto the floor (Figure 2C). Since the number of stabilized flights observed in the case of UB was very small (n = 3), no definite conclusions could be reached about the effects of this condition on the stabilization times, and these data were therefore removed from the statistical analysis. The lighting conditions significantly affected the stabilization times (p < 0.05; F = 5.5942): the mean stabilization time was approximately 50ms longer in the CB environment than that recorded in the two conditions with overhead lighting ( ∆Stab (ms): CB = 222.75;

T L (U T & CT ) = 169.0465).
A large number of stabilized flights occurred with overhead lighting and either textured or uniform lateral walls (Figure 2D). In the BL conditions, hoverflies produced poorer stabilization performances than in the TL conditions (p < 0.001; Chi 2 = 46.144), and they were almost unable to prevent themselves from continuing to fall in the uniform environment (UB). Hoverflies surrounded by textured lateral walls, i.e., in conditions CT and CB, achieved better performances than under the same lighting condition with uniform walls, i.e., in conditions UT and UB, respectively (p < 0.001; Chi 2 = 23.463). It was also observed that the stabilization rates were similar between Halogen-LED and control experiments (Effect of LED light: p = 0.11; Chi 2 = 2.5727; interaction with light position: p = 0.27; Chi 2 = 1.2127) whereas the effects of light position on stabilization time depended on the kind of light (Effect of LED light: p = 0.49; Chi 2 = 0.4730; interaction with light position: p < 0.01; Chi 2 = 7.52). It can be seen from figure 2C that this effect was mainly observed on the CT condition, where no differences seemed likely to occur. The differences were probably due to the fact that different populations were tested in the control and Halogen/Led experiments.

Discussion

In this study, it was attempted for the first time to investigate the impact of the light source position in hoverflies' visual field on their ability to stabilize their flight. In previous studies, it was suggested that the light gradient generated by an artificial horizon may impact insects' attitude perception processes via a mechanism called the Dorsal Light Response [START_REF] Hengstenberg | Visual Motion and its Role in the Stabilization of Gaze[END_REF][START_REF] Goulard | Behavioural evidence for a visual and proprioceptive control of head roll in hoverflies (episyrphus balteatus)[END_REF]. Using a free fall procedure, we reported that hoverflies starting to fly in an unsteady initial state were found to be able to recover stabilized flight efficiently only in situations where the light came from above. In addition to the crucial position of the light source, the optic flow information generated during a free fall may also be used by hoverflies to ultimately avoid crashing [START_REF] Hengstenberg | Mechanosensory control of compensatory head roll during flight in the blowfly Calliphora erythrocephala[END_REF][START_REF] Expert | Flying over uneven moving terrain based on optic-flow cues without any need for reference frames or accelerometers[END_REF][START_REF] Julien | Optic flow-based collision-free strategies: From insects to robots[END_REF][START_REF] Goulard | Modeling visual-based pitch, lift and speed control strategies in hoverflies[END_REF], but these cues probably do not suffice for stabilization purposes. The light gradient probably provides hoverflies with a means of estimating their absolute orientation in the environment in order to control their attitude, as found to occur in locusts [START_REF] Goodman | The Role Of Certain Optomotor Reactions In Regulating Stability In The Rolling Plane During Flight In The Desert Locust , Schistocerca Gregaria[END_REF]. This static cue providing a subjective vertical reference value about of the external world [START_REF] Hengstenberg | Localization and Orientation in Biology and Engineering[END_REF] would certainly require robust visual processing integrating the lighting information over the whole or most of the spherical field of view. To investigate in greater detail the extent to which the light gradient is actually involved, LED panels might be a useful means of finely controlling the homogeneity of the illumination generated and the amplitude of the light gradient, and generating dynamic changes in the lighting conditions.

One of the main hypotheses put forward in previous studies on flies' sensorimotor reflexes is that they may depend only on movement perception processes and compensatory reflexes [START_REF] Taylor | Sensory System and Flight Stability: What does Insects Measure and Why?[END_REF]. The inputs originating from both visual structures, the compound eyes and the ocelli [START_REF] Krapp | Estimation of self-motion by optic flow processing in single visual interneurons[END_REF][START_REF] Egelhaaf | Neural encoding of behaviourally relevant visual-motion information in the fly[END_REF][START_REF] Parsons | A motion-sensitive neurone responds to signals from the two visual systems of the blowfly, the compound eyes and ocelli[END_REF][START_REF] Parsons | Sensor Fusion in Identified Visual Interneurons[END_REF], and from the halteres [START_REF] Sandeman | Head Movements in Flies (Calliphora) Produced by Deflexion of the Halteres[END_REF][START_REF] Nalbach | The halteres of the blowfly Calliphora[END_REF][START_REF] Hengstenberg | Controlling the fly's gyroscopes[END_REF][START_REF] Dickinson | Haltere-mediated equilibrium reflexes of th fruit fly, drosophila melanogaster[END_REF], which are fused together non-linearly [START_REF] Huston | Nonlinear Integration of Visual and Haltere Inputs in Fly Neck Motor Neurons[END_REF],

may compensate for a large range of disturbances [START_REF] Schwyn | Interplay between Feedback and Feedforward Control in Fly Gaze Stabilization[END_REF]. A system of this kind is liable, however, to be subject to accumulated errors during flight, resulting in a drift in the attitude control process and eventually in crashing. The results obtained here therefore indicate that the DLR may play an important role by providing a reliable time-invariant vertical reference frame which may be used to complement the insects' motion-based reflexes. However, the initial position of hoverflies with their legs dangling may have decreased the ability of the chordotonal organs (organs acting as pressure sensors [START_REF] John | Mechanosensation and adaptive motor control in insects[END_REF] linked to postural reflex in insects [START_REF] Horn | Positional Head Reflexes and the Role of the Prosternal Organ in the Walking Fly, Calliphora erythrocephala[END_REF][START_REF] Kress | Head and body stabilization in blowflies walking on differently structured substrates[END_REF] which were

Figure 1 :

 1 Figure 1: Experimental setup.(A) The setup used in this study was first presented in [8]. In the present version, a white backlit LED panel was added to illuminate the box from below. (B) Four environments were tested: Contrasting textured lateral walls with Top lighting (CT), Contrasting textured lateral walls with Bottom lighting (CB), Uniform lateral walls with Top lighting (UT), and Uniform lateral walls with Bottom lighting (UB).

Four

  different conditions were tested (fig 1B): Contrasting textured lateral walls with Top lighting (CT), Contrasting textured lateral walls with Bottom lighting (CB), Uniform lateral walls with Top lighting (UT) and Uniform lateral walls with Bottom lighting (UB). In addition to these conditions, two control experiments were conducted in which the CT condition was compared with a CB condition with white halogen bottom lighting instead of the LED panel to check whether the LED lighting condition affected the hoverflies' performances. The texture on the walls consisted of a randomly generated chequerboard (20x20 squares 4cm 2

Figure 2 :

 2 Figure 2: (A) Boxplot of the wingbeat triggering times (ms). (B) Bar plot of the crash rates. (C) Boxplot of the stabilization times (ms). (D) Bar plot of the stabilized flight rates. Boxes are composed of 1st, 2nd and 3rd quartiles, and whiskers correspond to extreme data, amounting to no more than 1.5 times the interquartile distance. Significance code, p value: 0 < * * * < 0.001 < * * < 0.01 < * < 0.05

Funding

No funding support was involved in this research.

stimulated in this situation only by the legs´weight) to estimate their orientation with respect to the gravity experienced prior to the fall and we therefore cannot rule out the latter hypothesis.

In conclusion, this initial study shows that the position of the light source plays an important role in hoverflies' flight stabilization processes. The results presented here suggest that both static (light source, DLR) and movement (optic flow) cues are probably involved [START_REF] Hengstenberg | Localization and Orientation in Biology and Engineering[END_REF]. A further question which arises here is how these two visual processes (the DLR and OF-based processes) are fused together to ensure robust flight stabilization under natural conditions, as previously suggested in the case of optomotor responses [START_REF] Trischler | Chasing behaviour and optomotor following in free-flying male blowflies: flight performance and interactions of the underlying control systems[END_REF]. Future studies in which conflicting situations are generated would probably help to understand how these different sensory processes are combined in dipterans' brains.

Aknowledgments

We are most thankful to Julien Diperi for his contribution to building the experimental setup, to Marc Boyron for developing the electronics on which all the work presented in this paper was based, and to Jessica Blanc for correcting and improving the English manuscript.

Ethics

No ethical assessments for animal research purposes or permission to carry out fieldwork were required for this study.

Data accessibility

Data are provided in Excel tables in the supplementary files (Data total.xls).

Competing interests

We have no competing interests to declare Authors´contributions RG designed and conducted the experiments, analysed the data, carried out the statistical analyses and drafted the manuscript; AV conducted experiments and helped to draft the manuscript;

JLV helped with the statistical analyses and with the drafting of the manuscript. SV designed the experiment, coordinated the study and helped to draft the manuscript. All the authors gave their final approval for publication of the manuscript and agree to be held accountable for the work performed herein.


Florence Nicol
email: florence.nicol@enac.fr

Statistical Analysis of Aircraft Trajectories: a Functional Data Analysis Approach

Keywords: curve clustering, principal component analysis, functional statistics, air traffic management

In Functional Data Analysis, the underlying structure of a raw observation is functional and data are assumed to be sample paths from a single stochastic process. When data considered are functional in nature thus infinite-dimensional, like curves or images, the multivariate statistical procedures have to be generalized to the infinite-dimensional case. By approximating random functions by a finite number of random score vectors, the Principal Component Analysis approach appears as a dimension reduction technique and offers a visual tool to assess the dominant modes of variation, pattern of interest, clusters in the data and outlier detection. A functional statistics approach is applied to univariate and multivariate aircraft trajectories.

I. INTRODUCTION

In many fields of applied research and engineering, it is natural to work with data samples composed of curves. In air transportation, aircraft trajectories are basically smooth mappings from a bounded time interval to a state space. The dimension of the state space may considerably increase if Quick Access Recorders (QARs) provide a full bunch of flight parameters. Most of the time, aircraft trajectories are observed on a fine grid of time arguments that span the time interval. The size and the dimension of the observed samples are usually important, especially if the flight data recorders are used. Data collected in air transportation thus present some characteristics of big data: complexity, variety and volume. These characteristics are inherent to air traffic and require using specific statistical tools that take into account the diverse and complex nature of data and efficient numerical algorithms.

In Air Traffic Management (ATM), analyzing aircraft trajectories is an important challenge. A huge amount of data is continuously recorded (flight data recorder, maintenance softwares, Radar tracks) and may be used for improving flight, as well as airport safety. For instance, trajectories coming from flight data recorders might help the airlines to identify, measure and monitor the risk of accidents or to take preventive maintenance actions. On airports, landing tracks observations may indicate bad runway or taxiway conditions. Therefore, it is of crucial importance to propose relevant statistical tools for visualizing and clustering such kind of data, but also for exploring variability in aircraft trajectories.

Aircraft trajectories, that are basically mappings defined on a time interval, exhibit high local variability both in amplitude and in dynamics. Because of the huge amount of data, visualizing and analyzing such a sample of entangled trajectories may become difficult. A way of exploring variability is then to identify a small number of dominant modes of variation by adapting a Principal Component Analysis (PCA) approach to the functional framework. Some of these components can help to visualizing how major traffic flows are organized. This approach can also address the aircraft trajectories clustering that is a central question in the design of procedures at takeoff and landing. Moreover, identifying atypical trajectories may be of crucial importance in aviation safety. Resulting clusters and outliers may be eventually described relatively to other variables such as wind, temperature, route or aircraft type.

In this study, we will focus on Functional Principal Component Analysis (FPCA) which is a useful tool, providing common functional components explaining the structure of individual trajectories. First, in Section II and III, the state of the art and the general framework for functional data analysis are presented. Next, in Section IV, the PCA approach is generalized to the functional context. The registration problem is then considered when phase variation due to time lags and amplitude variation due to intensity differences are mixed. Finally, in Section V, FPCA is applied to aircraft trajectories that can be viewed as functional data.

II. PREVIOUS RELATED WORKS

Most of the time, aircraft trajectories are observed on a fine grid if time arguments that span the time interval. Data are first sampled then processed using multivariate statistics. While simple, this process will forget anything about the original functional dependency. Most of studies conducted on air traffic statistics make use of the sampled data only as is proposed in [START_REF] Eckstein | Data driven modeling for the simulation of converging runway operations[END_REF] and forget all about their functional nature, dropping some extremely valuable information in the process. One of the most salient shortcoming of the discrete samples methods is they do not take into account with the correlation in the data while functional data exhibit a high level of internal structure and intrinsic characteristics (geometry of trajectories). Moreover, as noted in [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF], standard methods of multivariate statistics have became inadequate, being plagued by the "curse of dimensionality". In a standard multivariate approach, a PCA is performed on matrix data in which the number of variables may be much more important than the number of individuals. As a result, statistical methods developed for multivariate analysis of random vectors are inoperative and trying to crudely apply traditional statistical algorithms on this kind of data may induce some severe numerical instabilities.

The quite recent field of functional statistics [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF] [3] provides a more adequate framework for dealing with such data that are assumed to be drawn from a continuous stochastic process taking its value in an Hilbert space. Data are no longer point values but the complete trajectories, all statistical procedure being performed on them. A major asset of working with functional data instead of points is the ease of adding a priori information by carefully selecting the Hilbert space. In air transportation, few studies using the functional framework have been carried out.

In [START_REF] Gregorutti | Grouped variable importance with random forests and application to multiple functional data analysis[END_REF], random forest for functional data are used for minimizing the risk of accidents and identifying explanatory factors in the context of aviation safety. This approach is not suitable to visualizing how major traffic flows are organized. In [START_REF] Puechmorel | Entropy Minimizing Curves with Application to Flight Path Design and Clustering[END_REF] [START_REF] Nicol | Unsupervised curves clustering by minimizing entropy: implementation and application to air traffic[END_REF], a new approach based on entropy minimization and Lie group modeling is presented, in which the geometry of trajectories are taken into account to cluster the traffic in groups of similar trajectories. Although this approach deals with the aircraft trajectories clustering, the objective is quite different. Indeed, this metod is intended to be a part of a future automated trajectory planner. Given a sample of planned trajectories, the classification algorithm creates clusters such that the mean line of each of them is similar to an airspace route. Geometrical constraints have then to be considered.

In [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF], a FPCA was performed on a sample of unidimensional aircraft trajectories, especially trajectory altitudes. This approach generalizes the standard multivariate principal component analysis described in [START_REF] Eckstein | Data driven modeling for the simulation of converging runway operations[END_REF] to the functional context. In the following, this approach is extended to the multivariate FPCA (MFPCA), in which we want to study the simultaneous modes of variation of more than one function. Particularly, the simultaneous statistical analysis of the longitude and latitude coordinates may give some insights on the nowadays traffic and then allow to forecast the expected one.

III. DEALING WITH RANDOM FUNCTIONS

A. Problem statement Functional Data analysis (FDA) deals with the study of infinite dimensional objects with a time or spatial structure to be processed, such as curves or images. This point of view differs from standard statistical approaches, the underlying structure of a raw observation being functional. Rather than on a sequence of individual points or finite-dimensional vectors as in a classical approach, we focus on problems raised by the analysis of a sample of functions. Functional data x 1 (t), . . . , x n (t) are the observations of a sample of n independent and identically distributed random functions X 1 (t), . . . , X n (t) that are assumed to be drawn from a continuous stochastic process X={X(t), t ∈ J}, where J is a compact interval. It makes sense to interpret functional data as n realizations of the stochastic process X, often assumed with values in a Hilbert space H, such as L 2 (J), the space of square integrable functions defined on the interval J. The associated inner product for such functions is x, y = x(t)y(t)dt and the most common type of norm, called L 2 -norm, is related to the above inner product through the relation x 2 = x, x . In a functional context, equivalence between norms fails and the choice of semi-metrics is driven by the shape of the functions, as noted in [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF]. For instance, semi-metrics based on derivatives suppose that the functions are not too rough.

Let X be a square integrable functional variable with values in the separable Hilbert space H. As noted in [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF], we can define a few standard functional characteristics of the random function X, such as the theoretical mean function and the theoretical covariance function, for s, t ∈ J,

µ(t) = E [ X(t) ] , (1) σ
(s, t) = E [ X(s)X(t) ] -E [ X(s) ] E [ X(t) ] , (2)
that play a crucial role in FPCA as we will see in Section IV.

In the following, we will assume that X is centered, that is µ = 0, otherwise, subsequent results refer to X -µ. From ( 1) and ( 2), we can derive the equivalent empirical characteristics. Note that no notion of probability density exists in the infinite dimensional Hilbert space as mentioned in [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF].

B. Trajectories smoothing

Usually, in practice, functional data, such as position and speed measurement, are observed discretely: we only observe a set of function values on a set of arguments that are not necessarily evenly space times or the same for all functions. Some preprocessing of the discretized data has to be made in order to recover the functional statistics setting, especially when observations are noisy. Most procedures developed in FDA are based on the use of interpolation or smoothing methods in order to estimate the functional data from noisy observations [START_REF] Ramsay | Functional Data Analysis[END_REF]. This problem can be solved by representing a trajectory as a linear combination of known basis function expansions such as a Fourier basis, wavelets or spline functions. Functional data are estimated by their projections onto a linear functional space spanned by K known basis functions ψ 1 , . . . , ψ K such as

x i (t) = K k=1 θ ik ψ k (t) = θ T i ψ(t), (3)
where the unknown coefficient vectors θ i = (θ i1 , . . . , θ iK ) T have to be estimated from the data and ψ(t) denotes the vectorvalued function (ψ 1 (t), . . . , ψ K (t)) T . Let us consider a set of sampled trajectories {(y ij , t ij ), i = 1, . . . , n, j = 1, . . . , N i } where y ij and t ij are the respective j-th sample position and time on the i-th trajectory. The argument values t ij may be the same for each recorded trajectory or also vary from one trajectory to another one. For simplicity, we will assume that the functional data are observed on the same time grid t 1 , . . . , t N , usually equally spaced. The expansion coefficient vector (θ i ) is the solution of the following least squares minimization problem

min θi j=1,...,N y ij -θ T i ψ(t j ) 2 = y i -Ψθ i 2 , (4)
where y i is the vector of the observed functional data and Ψ is the N × K matrix containing the values ψ k (t j ).

Note that this representation in a truncated basis functions takes into account the functional nature of the data and makes it possible to discretize the infinite dimensional problem by replacing the functional data x i (t) by its coefficient vector θ i , i = 1, . . . , n. While a probability density notion on an infinite dimensional Hilbert space cannot be defined [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF], the expansion of the curves on a truncated Hilbert basis allows to fit a distribution on the coefficient vectors. Usually, multivariate statistical procedures are next performed on the set of coefficients such as clustering techniques.

The choice of the number K of basis functions depends on the complexity of the curves. The larger is K in the expansion, the better is the fit but we then may capture undesirable noise. If K is too small, we may increase smoothness and some important characteristics of the functions may be vanished. Fixing the dimension of the model is not easy and a major drawback is due to the fact that the degree of smoothing is driven by the discrete choice of the parameter K. We can get better results by using roughly penalty approaches [START_REF] Ramsay | Functional Data Analysis[END_REF].

IV. A PRINCIPAL COMPONENT APPROACH Multivariate Principal Component Analysis is a powerful exploratory statistical method which synthesizes the quantity of data information by creating new descriptors in limited number [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF] [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]. FPCA was one of the first methods of multivariate analysis that has been generalized to a functional setting. As for the covariance matrix in the multivariate standard case, the covariance function of functional variables are difficult to interpret and FPCA goals to analyze the variability of the functional data around the mean function in an understandable manner. By approximating infinite-dimensional random functions by a finite number of random score vectors, FPCA appears as a dimension reduction technique just as in the multivariate case and cuts down the complexity of the data. For this reason, this approach is commonly used in FDA.

A. Generalization to the infinite-dimensional case

Let X 1 , . . . , X n be a sample of independent centered random functions. One wants to find weight functions γ i that preserve the major variation of the original sample. The criterion is then the sample variance of the projections of the random functions X 1 , . . . , X n into the weight functions, called principal component functions. These principal component functions are the solution of the maximizing problem:

max γi∈H 1 n n j=1 X j , γ i 2 , (5)
under the constraint:

γ i , γ k = δ ik , k ≤ i, i = 1, . . . , n. (6)
At each step, each principal component function represents the most important mode of variation in the random functions. The orthogonality constraint then provides an orthogonal basis for the linear subspace spanned by the random functions sample. The solutions are obtained by solving the Fredholm functional eigenequation that can be expressed by means of the sample covariance operator Γ induced by the sample covariance function σ:

Γ n v(t) = J σ n (s, t)v(s)ds (7) = 1 n n j=1 X j , v X j (t), v ∈ H. (8)
such that

Γγ i (s) = λ i γ i (s), s ∈ J. (9)
The principal component functions γ 1 , . . . , γ n are then the eigenfunctions of Γ n , ordered by the corresponding eigenvalues

λ 1 ≥ λ 2 ≥ • • • ≥ λ n ≥ 0. The projections A ij = γ i , X j , j = 1, .
. . , n are random variables, called principal component scores of X j into the γ i -direction [START_REF] Ramsay | Functional Data Analysis[END_REF]. These scores are centered, uncorrelated random variables accross j with variance equal to λ i . Another important property for FPCA involves the best L-term approximation property, meaning that the truncated expansion L i=1 A ij γ i is the best approximation of X j with a given number L of components in the sense of the mean integrated error. Because each functional variable X j admits the empirical Karhunen-Loève decomposition,

X j (t) = n i=1 A ij γ i (t), j = 1 . . . , n, (10)
the random scores A ij = γ i , X j can be interpreted as proportionality factors that represent strengths of the representation of each individual trajectory by the ith principal component function. Furthermore, FPCA provides eigenfunction estimates that can be interpreted as "modes of variation". These modes have a direct interpretation and are of interest in their own right. They offer a visual tool to assess the main directions in which functional data vary. As in the multivariate case, pairwise scatterplots of one score against another may reveal patterns of interest and clusters in the data. In addition, these plots may also be used to detect outliers and explain individual behavior relatively to modes of variation.

As in the multivariate PCA, we can easily measure the quality of the representation by means of the eigenvalue estimators. The ith eigenvalue estimator λ i measures the variation of the scores into the γ i -direction. The percentage of total variation τ i explained by the ith principal component and the cumulative ratio of variation τ C L explained by the first L principal components are then computed from the following ratio

τ i = λ i n i=1 λ i , τ C L = L k=1 λ k n i=1 λ i . ( 11
)
The amount of explained variation will decline on each step and we expect that a small number L of components will be sufficient to account for a large part of variation. Determining a reasonable number L of components is often a crucial issue in functional analysis. Indeed, choosing L = n components may be inadequate and high values of L are associated with high frequency components which represent the sampling noise. A simple and fast method to choose the dimension L is the scree plot that plots the cumulated proportion of variance explained by the first L components against the number of included components L. Alternative procedures to estimate an optimal dimension can be found in [START_REF] Kneip | Nonparametric estimation of common regressors for similar curve data[END_REF] and [START_REF] Besse | Pca stability and choice of dimensionality[END_REF].

B. Estimation

Several estimation methods of scores and principal component functions were developed for FPCA and asymptotic results was studied in [START_REF] Dauxois | Asymptotic theory for the principal component analysis of a vector random function: Some applications to statistical inference[END_REF]. The earliest method applied to discretized functional data to a fine grid of time arguments is based on numerical integration or quadrature rules [14] [15]. Numerical quadrature schemes can be used to involve a discrete approximation of the functional eigenequation ( 9)

Σ n W γ m = λ m γ m , (12)
where Σ n = ( σ n (t i , t j )) i,j=1,...,N is the sample covariance matrix evaluated at the quadrature points and W is a diagonal matrix with diagonal values being the quadrature weights. The solutions γ m = ( γ m (t 1 ), . . . , γ m (t N )) are the eigenvectors associated with the eigenvalues λ m of the matrix Σ n W . The eigenvectors γ m form an orthonormal system relatively to the metric defined by the weight matrix W . When the weight matrix W is not the identity matrix, an orthonormalization correction is needed using Gramm-Schmidt procedure. We can express the functional eigenequation in an equivalent symmetric eigenvalue problem

W 1/2 Σ n W 1/2 u m = λ m u m (13)
under the constraint:

u T l u m = δ lm , l, m = 1, . . . , N. (14)
where u m = W 1/2 γ m . Note that, if the discretization values t j are closely spaced, the choice of the interpolation method should not have a great effect compared to sampling errors, even if the observations are corrupted by noise [START_REF] Ramsay | Functional Data Analysis[END_REF].

A more sophisticated method is based on expansion of functional data on known basis functions such as a Fourier basis or spline functions as described in Section III. This method takes into account the functional nature of the data and makes it possible to discretize the problem by replacing the functional data x i (t) by its coefficient vector θ i , i = 1, . . . , n. The sample covariance function of the projected data

σ n (s, t) = 1 n n i=1 x i (s) x i (t) = ψ(s) T Θψ(t), (15)
can be expressed by means of the K × K matrix Θ =

where W = ( ψ i , ψ j ) i,j=1,...,K is the matrix of the inner products ψ i , ψ j = ψ i (t)ψ j (t)dt of the basis functions.

The solutions b m are then the eigenvectors associated with the eigenvalues λ m of the matrix ΘW . The orthonormality constraints on the principal components functions satisfy

b T l W b m = δ lm , l, m = 1, . . . , K. (17)
Note that this method looks like the discretization method for which the coefficient vectors θ i = (θ i1 , . . . , θ iK ) T play the role of the discretized functional data. FPCA is then equivalent to a standard multivariate PCA applied to the matrix of coefficients with the metric defined by the inner product matrix W = ( ψ i , ψ j ) i,j=1,...,K .

C. The registration problem

The process of registration, well known in the field of functional data analysis [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] [17] [START_REF] Ramsay | Functional Data Analysis[END_REF], is an important preliminary step before further statistical analysis. Indeed, a serious drawback must be considered when functions are shifted, owing to time lags or general differences in dynamics. Phase variation due to time lags and amplitude variation due to intensity differences are mixed and it may be hard to identify what is due to each kind of variation. This problem due to such mixed variations can hinder even the simplest analysis of trajectories. Firstly, standard statistical tools such as pointwise mean, variance and covariance functions, may not be appropriate. For example, a sample mean function may badly summarize sample functions in the sense that it does not accurately capture typical characteristics. In addition, a FPCA procedure applied to the unregistered curves will produce too many principal components, some of them being not of interest for the analysis of the variability of the curves. In addition, phase variation may influence the shape of the principal component functions that may not be representative of the structure of the curves. Finally, the scores may present a kind of correlation.

A registration method consists in aligning features of a sample of functions by non decreasing monotone transformations of time arguments, often called warping functions. These time transformations have to capture phase variation in the original functions and transform the different individual time scales into a common time interval for each function. Generally speaking, a non decreasing smooth mapping

h i : [a, b] → [c i , d i ], with [c i , d i ]
the original time domain of the trajectory, is used to map each trajectory y i to a reference trajectory x, usually called target or template function, already defined on [a, b]. In this way, remaining amplitude differences between registered (aligned) trajectories y i • h i can be analyzed by standard statistical methods. The choice of a template function is sometimes tricky and it may be simply selected among the sample trajectories as a reference with which we want to synchronize all other trajectories. Note that warping functions h i have to be invertible so that for the same sequence of events, time points on two different scales correspond to each other uniquely. Moreover, we require that these functions are smooth in the sense of being differentiable a certain number of times.

Most of literature deals with two kinds of registration methods: landmark registration and goodness-of-fit based registration methods. A classical procedure called marker or landmark registration aims to align curves by identifying locations t i1 , . . . , t iK of certain structural features, such as local minima, maxima or inflexion points, which can be found in each curve [START_REF] Bookstein | Morphometric Tools for Landmark Data: Geometry and Biology, ser. Geometry and Biology[END_REF] [19] [START_REF] Gasser | Searching for structure in curve sample[END_REF]. Curves are then aligned by transforming time in such a way that marker events may occur at the same time t 01 , . . . , t 0K , giving h i (t 0k ) = t ik , k = 1, . . . , K. Complete warping functions h i are then obtained by smooth monotonic interpolation. While this non-parametric method is able to estimate possibly non-linear warping functions, marker events may be missing in certain curves and feature location estimates can be hard to identify. Finally, phase variation may remain between too widely separated markers. An alternative method is based on goodness-of-fit by minimizing distance between registered trajectories and a template trajectory, with possible inclusion of a roughness penalty for h i [START_REF] Ramsay | Curve registration[END_REF] [START_REF] Ramsay | Estimating smooth monotone functions[END_REF]. Note that this latter registration method, as well as landmark registration are implemented in softwares R and Matlab [START_REF] Ramsay | Functional data analysis with R and Matlab[END_REF] and can be downloaded through the website [START_REF]Functional Data Analysis[END_REF].

V. APPLICATION TO AIRCRAFT TRAJECTORIES A. The aircraft trajectory dataset

We now apply the previously described FPCA technique to a 161 aircraft trajectory dataset. These data consist of radar tracks from Paris Charles de Gaulle (CDG) to Toulouse Blagnac airports recorded during two weeks. Most of the aircrafts are Airbus A319 (20%), A320 (18%) and A321 (33%), followed by Boeing B733 (15%), B463 (8%) a member of British Aerospace BAe 146 family and AT type (6%). Radar measurements are observed in the range of 4-6960 seconds at 4 seconds intervals. The assumption that all trajectories are sample paths from a single stochastic process defined on a time interval is clearly not satisfied in the case of aircrafts: departure times are different, even on the same origin-destination pair and the time to destination is related to the aircraft type and the wind experienced along the flight. Without loss of generality, we will assign a common starting time 0 to the first radar measurement of the flights. Trajectories in Figure 1 exhibit high local variability and may be studied by using a FPCA approach. As observed raw data were passed through preprocessing filters, we get radar measurements at a fine grid of time arguments with few noise. We have then used the discretization method described in Section IV.

B. Multivariate FPCA

We now apply the FPCA procedure to multidimensional trajectories. Each trajectory data f i (t) = (x i (t), y i (t)), i = 1, . . . , n, collected over time are effectively producing two dimensional functions over the observed intervals [0, T i ]. Trajectories have been registered by using the landmarks used in [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF] for the univariate altitude trajectories. Figure 2 displays the first four principal components for the latitude and longitude trajectories after the overall mean has been removed from each track. The first component in X and Y -coordinates explain 58.7% of total variation whose 98% is due to the longitude trajectories Y (t). We can visualize this effect on the overall mean function in Figure 3 second principal component is orthogonal to the first one, the corresponding mode of variation is less important and accounts for 14.7% of total variation. The contributions of both coordinates are of the same importance, with 48% and 52% of total variation respectively explained by X(t) and Y (t). In Figure 2, we can observe an overall effect due to the X(t) trajectories increasing with time and a distortion in the timing for the Y (t) trajectories. In Figure 3, we can visualize that the closer we get to Toulouse airport, the more aircraft trajectories are separated relatively to the X-coordinate. Moreover, the separation between the arrivals at Toulouse airport are slightly inflated relatively to the Y -coordinate. We call this effect the landing effect (PC2). The third component accounts for 12.9% and the main contribution comes from the X-coordinate with 86%. This component depicts an overall effect relatively to the X-coordinate that separates the two routes, immediately after the take-off from Paris CDG airport. We call this effect the separation effect (PC3). Finally, the fourth principal component accounting for 6% of the total variation, whose 66% is explained by the X-coordinate, highlights an inversion of route, probably due to a change of take-off procedures at Paris CDG airport or landing procedures at Toulouse airport. We call this effect the change effect (PC4). A k-means clustering is next performed on the score matrix. In Figure 4, we can visualize the mean cluster trajectories for three and five clusters. The first cluster (blue line) contains all aircraft types except the AT type while the third one (red line) is mainly composed of AT type. The mean trajectory of the first cluster displays the overall flight paths from Paris VI. CONCLUSION AND FUTURE WORKS FPCA is a powerful tool to analyze and visualize the main directions in which trajectories vary. We have successfully applied this technique to analyze aircraft trajectories and it can be easily extended to the multivariate case. FPCA has many advantages. By characterizing individual trajectories through an empirical Karhunen-Loève decomposition, FPCA can be used as a dimension reduction technique. Moreover, rather than studying infinite-dimensional functional data, we can focus on a finite-dimensional vector of random scores that can be used into further statistical analysis such as cluster analysis.

The FPCA approach seems promising, as indicated by the results obtained on a real data set. However, the registration problem remains crucial because the assumption that all trajectories are sample paths from a single stochastic process is not satisfied and may be complex in the case of multidimensional aircraft trajectories. In this work, we have used a landmark registration technique. In future works, we will use more sophisticated procedures such as arclength parametrization.

Moreover, we should add heading and velocity information by combining functional data and vector of data, inducing an extra level of complexity.

  covariance matrix of the coefficient vectors. Consider now the basis expansion of the eigenfunctions γ m (s) = b T m ψ(s) where b m = (b m1 , . . . , b mK ) T is the unknown coefficient vector to be determined. This yields the discretized eigenequation ΘW b m = λ m b m ,

Figure 1 .

 1 Figure 1. Trajectories from Paris CDG airport to Toulouse airport.

Figure 2 .

 2 Figure 2. The first four principal component functions for the latitude trajectories X(t) and the longitude trajectories Y (t).

Figure 3 .

 3 Figure 3. The effects on the mean aircraft trajectory (solid curve) of adding (red curves) and substracting (blue curves) a multiple of each of the first four principal components.

Figure 4 .

 4 Figure 4. Mean cluster trajectories and the overall mean (black curve).

  by adding and subtracting a suitable multiple of the first principal component for each coordinate. This component quantifies an overall decrease in longitude that we can call overall effect (PC1) between the two different routes from Paris (CDG) to Toulouse airports, more and more important when one moves towards Toulouse airport. Aircrafts with high negative scores would show especially above-average tracks, mainly due to the Y -coordinate. As the

				PC1 58.7 %					PC2 14.7 %
		100						100
		50						50
		0						0
	Y(t) (en Nm)	-50					Y(t) (en Nm)	-50
		-100						-100
		-150						-150
		-200						-200
		50	60	70	80	90	100	50	60	70	80	90	100
				X(t) (en Nm)					X(t) (en Nm)
				PC3 12.9 %					PC4 6 %
		100						100
		50						50
	Y(t) (en Nm)	-50 0					Y(t) (en Nm)	0 -50
		-100						-100
		-150						-150
		-200						-200
		50	60	70	80	90	100	50	60	70	80	90	100
				X(t) (en Nm)					X(t) (en Nm)

TABLE I .

 I CONTINGENCY TABLE OF THE COUNTS

	Aircraft type	Cluster 1	Cluster 2	Cluster 3
	A319	15	18	0
	A320	14	14	1
	A321	25	28	0
	AT	2	0	8
	B463	10	0	2
	B733	22	1	2

TABLE II .

 II CONTINGENCY TABLE OF THE COUNTS

	Aircraft type	Cluster 1	Cluster 2	Cluster 3	Cluster 4-5
	A319	9	16	0	8
	A320	10	13	0	6
	A321	18	13	0	6
	AT	0	0	8	2
	B463	10	0	2	0
	B733	17	0	1	6


Stéphane Florence Nicol

Curves Similarity Based on Higher Order Derivatives

come   L'archive ouverte pluridisciplinaire

Curves Similarity Based on Higher Order Derivatives

Florence Nicol †, Stephane Puechmorel ‡ Ecole Nationale de l'Aviation Civile Email: †nicol@recherche.enac.fr, ‡stephane.puechmore@enac.fr Abstract-In many applications, data originate from the observation of a phenomenon depending on time. Trajectories of mobiles fall within this category and receive an increasing attention as many connected objects have the ability to broadcast their positions. When the raw location is the value of interest, several statistical procedures exist to deal with analysis of trajectories. Depending on whether the geometrical shape or the time to position relation is relevant, one will use a parametrization invariant distance or a simple L 2 metric to assess the similarity between any two trajectories. However, it is sometimes advisable to use higher order information like velocity or acceleration, while retaining some kind of geometrical invariance. The purpose of the present work is to introduce a framework especially adapted to such a situation.

keywords-bundle metric, curve manifold, shape space.

I. INTRODUCTION

In many applications, the data of interest are measured values through time of a system in evolution. It is often the result of the observation of a physical phenomenon, obeying an underlying dynamics that may be unknown. As an extension, images can be modeled pretty much the same way, using two coordinates instead of a single time axis. Most algorithms designed to deal with such data rely on a sampled representation that is amenable to multivariate statistics.

Another approach was taken in the functional statistics framework, where the basic objects are mappings from a time interval to a given state space. Unfortunately, very little is known about probabilities in infinite dimension spaces, and one has to revert to finite dimensional representations during the implementation phase.

Several works were dedicated to the extension of classical multivariate algorithms to sample paths of Hilbert processes. As a starting point, data is first expanded on a truncated Hilbert basis [START_REF] Ramsay | Functional Data Analysis[END_REF], then, the vectors of expansion coefficients enter a standard finite dimensional analysis. A clever choice of the representation space and basis allows to take into account the prior knowledge about the studied process. Unfortunately, the dimension of the samples produced that way may be high, and varies with the geometric features of the sample paths. In particular, the presence of high curvature values will increase the number of expansion coefficients needed to keep a good approximation of the original function. In [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF], an expectation maximization (EM) functional clustering algorithm is presented with an adaptive basis in each group, yielding an efficient numerical method to deal with this issue.

Another class of methods relies on a non-parametric approach [3][4]. A recent work [START_REF] Boullé | Nonparametric Hierarchical Clustering of Functional Data[END_REF] pertaining to this approach presents a hierarchical clustering principle, with application to electric power consumption.

Finally, some algorithms use a shape manifold [START_REF] Mumford | The geometry and curvature of shape spaces[END_REF] representation in order to derive a metric between sample paths. This last class of methods has distinguished benefits, like the ability to focus only on the shape of the curves and forget about a specific parametrization. The major drawback is a high computational cost, that may preclude its use on large data sets. While basically designed for using the first derivatives, it is possible to consider higher order information, although the notion of parametrization invariance becomes less intuitive.

The purpose of the present work is to introduce a framework in which the higher order derivatives are explicitly taken into account, but with a well controlled notion of invariance. It was motivated by an application in civil aviation, that is the assessment of runway adherence using only radar tracks of landing and taxiing aircraft. In this context, a full parametrization invariance is not advisable, as it will remove an important part of the relevant information. On the other side, a raw Sobolev distance between curves will be fooled by the diversity of aircraft and will induce false alarms.

In section (II), a general approach to the question of curve similarity is addressed. The main idea is to split between the observation and the geometrical object on which it lies, yielding a metric that gather in a controlled way the contribution of the underlying shape and the extra information borne by the measured data. Starting with a model of curves based on differential geometry, the data can be thought as a section of a vector bundle with base space the curve. Using a riemanian metric on it will be the key to obtain a measure of similarity between any two samples, in the spirit of the works dedicated to distances between shapes [START_REF] Mumford | The geometry and curvature of shape spaces[END_REF].

In section (III), the application of the general framework to the case of landing aircraft tracks will be discussed. Finally, a conclusion will be drawn, with a view towards future work on application on real and simulated data prior to an operational use.

II. DISTANCE BETWEEN BUNDLES

A. Problem statement

The purpose of this section is to introduce a suitable state space for representing data that has both geometrical and cinematic features. It was motivated by the application that will be described later, where tracks of landing and taxiing aircraft must be clustered into homogeneous groups distinguished by the adherence of the runway. An example of a nominal (lower) and an abnormal (upper) track is given in Figure 1. The shape of the upper track is clearly different from the nominal one, and will not belong to the same cluster if one use a algorithm based for example on curvature. However, the divergence from the nominal curve may be due to low adherence condition or to a late turn. Only the first case must trigger a corrective action, that for the present situation is quite constraining: the runway has to be closed for the duration needed to perform an on-site adherence measurement. Using tangential acceleration will add a very discriminating feature since for the same shape, a high value means a high braking force, thus rejecting the hypothesis of low adherence. However, comparing raw velocity and acceleration induces two new issues:

• The aircraft type and airline procedures are varying, so that a shape comparison must be performed to minimize this effect; • The time span of the trajectories is not the same. A registration procedure is needed, and it is a quite difficult problem to solve. To summarize, neither parametrization invariant nor time dependent distances are fully satisfying. It is thus advisable to split the similarity computation into a purely geometrical part and a remainder that is not explained by shape variations.

All curves will be assumed to be smooth, that is indefinitely differentiable. This is not a real constraint in practical applications.

B. Curves as manifolds

In almost all applications, a curve is understood as a mapping γ from a real interval [a, b] to a state space, generally R p . Furthermore, only curves with nowhere vanishing first derivative will be considered to avoid possible singularities. In practice, this condition is never an issue. This is the standard framework in which functional data statistics takes place, and is well suited to problems where the information is contained in the mapping. As an example, if one wants to analyze the delays occurring in road or air traffic, the time to position mapping is the most relevant data. However, it is quite common to encounter cases where the shape of the image γ([a, b]) holds the discriminating features. It obviously the case in image recognition algorithms, but also in spectrometric data [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF], in biometric measurements (electroencephalogram, electrocardiogram) and generally speaking in all areas where the information will not change if the parametrization of the curves is changed. In the sequel, such a situation will be referred to as geometric data analysis (GDA).

The easiest approach to GDA is to let all curves be parametrized by arclength, which is defined as:

s : t ∈ [a, b] → t a γ (t) dt
The arclength has domain [0, l(γ)] with l(γ) the length of the curve. It comes at once that:

ds dt = γ (t)
so that taking the derivative with respect to s of a function of t can be done easily using the formula:

D s = 1 γ (t) D t
The arclength is related only to the shape of the image of γ.

In fact, let u ∈ [c, d] such that t = φ(u) with φ a strictly increasing smooth diffeomorphism from [c, d] to [a, b], it comes: s(u) = u c D u γ(t(u)) du = u c γ (t(u)) dt du du = t a γ (t) dt (1)
Many important features of the curve are naturally expressed with arclength: D s γ(s) yields the unit tangent vector T γ (s) while D ss γ(s) is the curvature at s. For curve similarity computations, the drawback of the arclength is that its domain depends on the length of the curve and varies with the curves. A scaled version η with domain [0, 1] is more convenient:

η = s/l(γ)
An obvious benefit of using η as a reference parametrization is that it solves the so-called registration problem, that consists of finding a common domain for all the functional samples [START_REF] Sangalli | Functional clustering and alignment methods with applications[END_REF]. It worth notice that in a sampled context, it is equivalent to have evenly spaced landmarks on the image of γ, as in [START_REF] Dryden | Statistical Shape Analysis, ser. Wiley Series in Probability & Statistics[END_REF].

It is worth noticing that the geometry of smooth curves with values in R p is entirely defined by the so-called Frenet frame (u 1 , . . . u p ) and its associated curvatures (κ 1 , . . . , κ p-1 ). For the sake of completeness, the procedure for finding them is given by

u 1 (t) = γ (t)/ γ (t) (2) ũi (t) = γ (i) (t) - i-1 j=1 γ (i) (t), u j (t) u j (t), i = 2 . . . p (3) u i (t) = ũi (t)/ ũi , i = 2 . . . p (4)
κ i (t) = u i (t), u i+1 (t) γ (t) , i = 1 . . . p -1 (5)
It is clear from the construction that changing the curve parametrization will keep the Frenet frame and the curvatures invariant: all information related to velocity, tangential acceleration and so on will be lost. In many applications, this is a major issue. From a mathematical standpoint, one possible geometric model for a curve is a one dimensional riemanian manifold. It will be assumed in the sequel that the reader is familiar with the basic concepts of differential geometry that may be found in any textbook on the subject [START_REF] Flaherty | Riemannian Geometry, ser. Mathematics: Theory & Applications[END_REF]. For shape analysis and recognition, one deals almost always with closed curves, that comply with the usual notion of manifold. However, when dealing with paths with distinct endpoints, the right model is a manifold with boundary. In the sequel, all curves will be represented that way.

Any real interval [a, b] has the structure of a trivial one dimensional manifold with boundary {a, b}. Tangent vectors are couples (t, u) where t is the basepoint in [a, b] and u is a one dimensional vector, that can be represented as a real number. One can think of a tangent vector (t, u) as the velocity at t of a point moving along the interval [a, b]. Boundary conditions impose u(a) > 0, u(b) < 0 . A vector field defined on it is just a smooth mapping u : t ∈ [a, b] → R such that all derivatives admit a limit to the right (resp. to the left) at a (resp. b). These limits will define the field at the boundary, and are not required to comply with the conditions satisfied by tangent vectors.

A curve γ : [a, b] → R p with nowhere vanishing derivative is an immersion from the manifold ]a, b[ to R p that can be extended to [a, b]. As such, it inherits a metric from the local embedding in R p by letting, for any t ∈]a, b[ and real numbers u, v, interpreted as tangent vectors:

g t (u, v) = γ (t) 2 uv
The Levi-Civita connection on ]a, b[ is given by:

∇ ∂t u = ∂ t u + 1 γ (t) 2 γ (t), γ (t) u
The first term corresponds to the intrinsic variation with respect to the parameter t, while the second is the part of the tangential acceleration coming from the geometry of the immersion. When the immersion parameter is chosen to be the arclength, the second term in the right hand vanishes as D s γ is the unit tangent vector at s and D ss γ is orthogonal to it. The Levi-Civita connection reduces thus to derivative with respect to arclength. The same applies for the scaled arclength η.

While curvature is a very important information for curves, it is a characteristic of the immersion and not an intrinsic feature of the trivial manifold [a, b]. It can be recovered from the immersed normal bundle N , whose elements are couples (t, v) with t ∈]a, b[ and v ⊥ γ (t). It is a vector bundle with base manifold [a, b] and typical fiber R p-1 , and its sections can be easily recovered using the Frenet frame u 1 (t), . . . , u p (t) introduced earlier, as any vector normal to γ (t) lies in span(u 2 (t), . . . , u p (t). A section of N is then a smooth mapping s : [a, b] → R p-1 , with s(t) the coordinates on the frame (u 2 (t), . . . , u p (t)).

The immersed normal bundle is the prototype of objects bearing vector information along a curve, and fulfilling the requirement of partial geometrical invariance mentioned earlier. For a single curve γ : [a, b] → R p , one can attach a vector v(θ) for each θ ∈ [a, b] that is interpreted as a sample at position γ(θ). Going back to the case of landing tracks, one can think of this vector information as the couple velocity/acceleration v(θ) = (u(γ(θ)), D t u(γ(θ)). It is important to note that the shape parameter θ is not related to time t, that is used to compute velocity and acceleration: at a given θ, u(θ), D t (θ) are the respective velocity and acceleration sampled at position γ(θ) on the trajectory. Different aircraft following the same curve but with different braking profiles will have different samples v(θ). Within the immersed bundle, they will be described by different sections. However, the geometric object underlying the sections, that is the base immersion γ will remain the same. When we let it vary, it appears that data can be compared at two well defined levels: the first one is the geometry and arises from the difference between the base immersions and the second is related to the sections themselves.

In a general setting, a immersed vector bundle will be defined as a vector bundle with base manifold [a, b] , typical fiber R n and whose sections are of the form s(t) = v(γ(t)), t ∈ [a, b], with γ : [a, b] → R p an immersion. As above, the geometry is encoded in γ, while the non-geometric information is described by the section.

The ability to deal with similarity between such objects relies on a notion of distance between them, that will be now introduced.

C. Distance between bundle sections

In order to simplify the derivations, all curves will be assumed to be immersions from [0, 1] to R p , using the η parameter. The derivation of a distance between immersed bundles sections will closely follow the principle underlying the construction of geometric distances between curves, as presented in [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF]. Let E 0 , E 1 be immersed vector bundles on respective immersions γ 0 , γ 1 with values in R p and with typical fibers R n . Let s 0 , s 1 be sections respectively of E 0 , E 1 , that will represent the vector samples along the respective curves γ 0 , γ 1 . An immersed path between s 0 and s 1 is a smooth mapping φ :

[0, 1] × [0, 1] → R p × R n such that: • For all s ∈ [0, 1], the mapping t ∈ [0, 1] → φ(s, t) is a smooth section the trivial bundle R p × R n → π R p ; • For all s ∈ [0, 1], π • φ(s, •) is a smooth immersion in R p ; • For all t ∈ [0, 1], φ(0, t) = s 0 (t), φ(1, t) = s 1 (t).
For a given s ∈ [0, 1], the mapping t ∈ [0, 1] → π • φ(s, t) defines an immersion from [0, 1] → R p that interpolates between γ 1 , γ 2 in the sense of [10][6]. This immersion defines in turn an immersed bundle, with typical fiber R n . It will be referred to as E s in the sequel. Finally, if a metric g s is available on each of the member of the family E t , it is possible to compute for a given φ a path length:

l(φ) = 1 0 1 0 g s ∂φ ∂s , ∂φ ∂s 1/2 dsdt (6)
Having this measure at hand, the distance between s 0 and s 1 is defined the infimum of the values l(φ) over all the admissible paths φ.

For clustering applications, it may be convenient to use the energy of a path φ that is defined to be:

E(φ) = 1 0 1 0 g s ∂φ ∂s , ∂φ ∂s dsdt
Paths minimizing the energy are the same as those minimizing the length. Since it saves a square root in the computation, the minimization is easier. However, if a distance is really needed at the end, it is enough to compute l(φ) on the minimizing path obtained. The way the interpolating bundles E s can be constructed will be deferred to a future work; however, normal bundles (or closely related ones) are quite natural in applications. It was the choice made for the classification of landing aircraft tracks.

A second question is the choice of the metrics g s . While out of the scope of the present paper, a requirement is that the family be smooth and natural in the sense of [START_REF] Kowalski | Natural transformations of riemanian metrics on manifolds to metrics on tangent bundles[END_REF]. For the application, an ad-hoc metric will be derived in the next section.

III. A WORKED EXAMPLE: SKID DETECTION

A. Problem statement

The problem of early detection of runway bad adherence has a great importance in airport operational management. When the runway or the taxiways are in bad condition, the only reliable procedure used nowadays is a direct measurement using an especially designed vehicle. Unfortunately, the runway has to be closed for the duration of the operation, which has a high cost both from the economic and traffic management point of view. Attempts where made to infer adherence from clustering of landing trajectories obtained by the surface surveillance means (radars), but due to the diversity of aircraft and airline procedures, a registration must be applied to curves, which is quite awkward to design. Since both the geometry of the landing trajectories and the deceleration law are important to make the right decision, the above theoretical framework seems to be ideally suited. As the phenomenon of interest stems from contact mechanics, it is worth starting with the underlying physics to derive a suited bundle metric.

B. From contact mechanics to curve similarity

Landing aircraft may experience slip during deceleration phase when the runway is in degraded conditions. It may result from icing, snow, bad runway surface state but also from pilot's actions, namely a too strong braking action or sharp turn. In this last case, it is not related to runway condition and must not trigger a maintenance action from the airport services.

Slip can be detected on-board by comparing wheel rotation rate with aircraft velocity and computing the so-called wheel slip factor:

λ = ω w -ω a max(ω w , ω a ) (7)
where ω w is the wheel angular velocity and ω a = V a /R w is the expected angular velocity that can be computed as the ratio of the aircraft velocity to the wheel radius. Please note that on the real vehicle, several wheels are used, and the λ coefficient has to be understood as a mean value. Furthermore, due to tire elasticity, λ is not zero even if there is no actual slip: this is due to the fact that when a traction or a braking force is applied, the rubber will stretch, resulting in the tire outer part actually traveling more or less than expected from rigid body dynamics. This information is not yet downlinked in real time to ground centers and thus cannot be used in the intended application. From the ground standpoint, λ cannot be observed without on-board information, but some aspects of the landing or taxiing aircraft behavior may still be inferred. It is assumed in the sequel that Coulomb's law for friction [START_REF] Popov | Contact Mechanics and Friction: Physical Principles and Applications[END_REF] is applicable, so that the contact force F c depends only on aircraft weight and tire/runway conditions:

F c ≤ µgM (8)
with M the aircraft mass, g the gravity of the Earth and µ the adhesion coefficient. Without slip, µ is equal to the static friction coefficient µ s and F c can be increased until it reaches the upper bound in [START_REF] Dryden | Statistical Shape Analysis, ser. Wiley Series in Probability & Statistics[END_REF]. At that point, slip occurs and µ drops to the value of the dynamic friction coefficient µ d . F c remains constant until it falls below µ d gM . In real world experiments, this simple behavior is no longer valid and one has to express µ as a function of λ, which can be found in [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]. Within this frame, the expression of the contact force is F c = µ(λ)gM , which is valid for both non-slip and slip case. Furthermore, in the case of aircraft, aerodynamics forces are exerted, with a net result of a braking force F a that adds to the actual brakes action, but does not contribute to the friction analysis. Putting things together, the equation of motion along the aircraft trajectory γ can be expressed as:

γ(t) = F a (t) M + µ(λ(t))g u ( 9
)
where u is a unit vector in the direction of the contact force F c . Without making additional assumptions, it is not possible to use [START_REF] Flaherty | Riemannian Geometry, ser. Mathematics: Theory & Applications[END_REF] for slip detection. However, if actions taken are assumed to be optimal, then F a and u will be collinear so as to maximize the net braking effect. The expression of the aircraft dynamics becomes:

γ(t) = (K(t) + µ(λ(t)g) u (10)
where the coefficient K(t) accounts for the aerodynamic braking force intensity. As aircraft must loose speed fast, µ will be close to the maximum at least during the landing and the beginning of taxi. The same applies for K, as it will not impair adherence. It can then be deduced that aircraft will try to keep the ratio between longitudinal and normal acceleration as high as possible. An observable measurement of slip condition will then be given by:

θ = arctan κ D t γ 3 D tt γ, D t γ ( 11
)
where γ is the aircraft trajectory and κ its curvature. It can further reduced to:

θ = arctan det (D t γ, D tt γ) D tt γ, D t γ (12)
using the well known property κ = det (D t γ, D tt γ) / D t γ 3 . In good runway conditions, the longitudinal acceleration will be high and nearly constant, at least in the first part of the landing trajectory. As a consequence, one can expect θ to be relatively small and be proportional to det (D t γ, D tt γ). Reciprocally, under slip conditions, a trade off has to be made between path following and deceleration: the angle θ will thus increase towards the limiting value ±π/2.

From the above discussion, it appears that θ makes sense as a weighting factor for curve comparisons.

C. An adapted metric

In the sequel, the symbol D t will stand for the partial derivative with respect to variable t. Higher order derivatives with respect to variables t 1 . . . t N will be written similarly as D t1t2...t N . Please note that the same variable may occur several times.

A smooth planar curve γ : [0, 1] → R 2 will be an immersion when the derivative D s γ is everywhere non vanishing in ]0, 1[. The set of such curves will be denoted by Imm([0, 1], R 2 ). Taking γ ∈ Imm([0, 1], R 2 ), its length is:

l(γ) = 1 0 D t γ(t) dt (13) It is invariant under parametrization change γ → γ • φ with φ : [0, 1] → [0, 1] a smooth diffeomorphism. Given a path Φ : [-, ] → Imm([0, 1], R 2
) that can be seen as a smooth mapping Φ : [-, ] × [0, 1] → R 2 , the variation of l(Φ(0, •)) can be computed :

D s | s=0 l(Φ(s, •)) = D s Φ(0, 1), T (1) -D s Φ(0, 0), T (0)
+ 1 0 D s Φ(0, t), κ(t) D t Φ(0, t) 2 N (t) dt (14)
with T (t), N (t) the respective unit tangent and normal vectors to the curve t → Φ(0, t) and κ(t) its unsigned curvature. The extension to more general immersions is quite straightforward [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF]. In the same reference, the variation formula ( 14) is used to derive a riemanian metric on the quotient space

Imm(S 1 , R 2 )/Diff ([0, 1], R 2 ).
In the present work, a similar approach will be taken. However, due to the fact that the slip condition must come into play as a weighting factor, it is not possible to keep invariance under change of parametrization. Furthermore, curves with vanishing second derivative must be excluded since the slip angle θ 11 is not defined at points where D t tγ(t) = 0. The last condition boils down to the requirement that the curve t ∈ [0, 1] → (γ, D t γ) be an immersion. The space of such objects will be denoted by Imm([0, 1], R 4 ).

As an arc tangent appears in the definition of θ, it is more convenient to use instead sin(θ) that has a simpler expression but otherwise similar behavior:

sin(θ(t)) = κ(t) D t γ(t) 2 D tt γ(t) = det (D t γ(t), D tt γ(s)) D t γ(t) D tt γ (16)
Definition 1. Let γ be a smooth curve. An admissible variation of γ is a smooth mapping

Φ : ] -, [→ R 2 , > 0, such that Φ(0, •) = γ(•) and ∀s ∈] -, [, Φ(s, 0) = γ(0), Φ(s, 1) = γ(1).
An admissible variation defines a tangent vector at γ: it is the smooth vector field t ∈ [0, 1] → D s Φ(0, t).

The expression (16) has a nice variational interpretation as indicated in the next lemma.

Lemma 1. Let γ : [0, 1] → R 2 be a smooth path and Φ an admissible variation of it. Let φ be a smooth path such that φ(0) = γ(1), φ(1) = γ(0). Then:

D s A(0) = 1 0 det (D s Φ(0, t), D t (0, t)) dt (17)
where A(s) is the net area enclosed by the loop Φ(s, •), φ for s ∈] -, [.

Using lemma 1, the integral : Turning back to bundles, if γ : [0, 1] → R 2 is the immersion describing the geometry of the problem (with the η parametrization) and v(η), v (η) the respective velocity and acceleration at position γ(η), the above metric can be adapted to yield a bundle metric. Without going into derivation detail, it can be expressed as: Looking at the expression (19) reveals a sum of a geometric distance between immersions in the sense of Mumford-Michor and a proper vector variation. This bundle metric is injected in the procedure for computing distances between immersed bundles sections [START_REF] Mumford | The geometry and curvature of shape spaces[END_REF] to yield the desired similarity measure that discriminates skid conditions.

IV. CONCLUSION AND FUTURE WORK

While quite developed from the theoretical standpoint, the work is still ongoing to assess performance in operational environment on real data. Unfortunately, there is no access to data correlated to adherence condition, as this information can only be obtained from direct measurements and is not communicated by airports authorities. The option taken was to develop a realistic taxi and landing simulator, that has been recently completed and will be released in open source. With the help of this tool, slipping and non-slipping trajectories can be simulated and the performance of the classification procedure assessed. On available landing data, and using an upper bound estimate of the distance based on a linear homotopy as an admissible path, promising results have been obtained. However, they cannot be matched against adherence conditions due to the aforementioned data availability issue. It is nevertheless expected, based on this experiment, that even the simplest linear homotopy procedure will outperform all the state-of-the-art algorithms investigated so far.

Figure 1 .

 1 Figure 1. Runway clearing trajectories

1 0

 1 |det (D t γ(t), D tt γ(s))| D t γ(t) D tt γ D t γ(t) dt(18)may be interpreted as the total infinitesimal area swept by the curve γ when moved in the direction D tt γ. This quantity is global indication of the slipping experienced along the path γ.

  g((u(η), u (η)), (v(η), v (η)) = u(η) N , v(η) N (1+κ 2 (η)) + det (D η γ(t), D ηη γ(s)) (19)