ENAC Florence Nicol

Baptiste Gregorutti

Safety Line

Cindie Andrieu

Curve clustering based on second order information: application to bad runway condition detection

Keywords: curve clustering, outlier detection, similarity measure, shape manifold, functional data, air traffic management, airport safety

In air transportation, a huge amount of data is continuously recorded such as radar tracks that may be used for improving flight as well as airport safety. However, all known statistical algorithms such as clustering procedures, even those based on functional data, are unable to distinguish between a safety critical flight and another one departing from standard behavior, but otherwise safe. It is the case in airport safety when radar measurements are used for detecting incidents on airport surface. In this paper, we propose a change of paradigm by switching from a functional data framework to a geometrical one by representing curves as points in a shape manifold. In this way, any intrinsic structure of the data that is amenable to geometry can be directly encoded in the representation space. Based on an extension of a classical distance between shapes, a new one is defined, that explicitly takes into account the second derivative and can be related to slippery. Its properties are investigated in a first part, then some results on datasets of synthetic and real trajectories are presented.

I. INTRODUCTION

I N the context of aviation safety, a crucial issue consists in assessing runway adherence condition and detecting incidents on airport surface by observing only the radar tracks of landing aircraft, e.g. trajectory deviation and abnormal deceleration may be good indicators of runway bad state. Indeed, landing aircraft must brake in a quite short time, putting a stringent condition on the adherence coefficient of the runways. Bad weather conditions like rain, snow and icing can dramatically lower it, increasing the landing distance and making maneuvers more difficult to perform. Rubber deposited by the wheels during braking may also impair the friction coefficient of the runway making it more slippery.

To estimate the runway adherence condition, a direct measurement of it on the pavement is usually performed: it implies sending a vehicle with a dedicated tool, which will interfere with the ongoing traffic. With the expected increase in air traffic, this procedure is particularly restrictive: based on recent studies, European air traffic was expected to grow, yielding an increase of 1.2 million flights in 2023, 14% more than in 2016 [START_REF]Annual Analysis of the EU Air Transport Market[END_REF]. In the context of the COVID-19 pandemic, the massive drop in air traffic and the traffic evolution significantly have reduced this problem. However, it still arises at hub airports and this procedure has a significant operational cost and environmental impact. Major framework programs such as SESAR (Single European Sky Air traffic management Research) in Europe aim to improve air traffic management (ATM) by investigating innovative solutions on management and traffic flow analysis. In particular, for assessing runway adherence condition, we need for reliable, continuous and consistent tools that not requiring runway closure. The objective is to develop spatio-temporal similarity measures for trajectories that will be used in clustering or outliers detection methods.

A. Motivations for Functional Data Analysis

The mathematical problem arising from the application presented above falls within the general framework of functional data statistics. Indeed, aircraft trajectories are functional objects mapping time to position, even if we will most of the time observe discretized samples of trajectories, such as radar measurements.

In the functional data framework, the commonly used algorithms developed for multivariate analysis are inoperative and may induce some numerical instabilities. Several works were dedicated to the extension of 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 enters a standard finite dimensional analysis. A clever choice of the representation space and basis allows to take into account the a priori knowledge about the studied process. In [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF], an EM functional clustering algorithm is presented with 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 [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice, ser. Springer Series in Statistics[END_REF], [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. 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.

In ATM, some works have already used functional data statistics, such as functional principal component analysis (FPCA), functional regression and random forest. FPCA is performed on aircraft trajectories in [START_REF] Nicol | Functional principal component analysis of aircraft trajectories[END_REF], [START_REF]Statistical Analysis of Aircraft Trajectories: a Functional Data Analysis Approach[END_REF], yielding some insights on the nowadays traffic and then allow to forecast the expected one. Such a functional decomposition is used in [START_REF] Hurter | Functional Decomposition for Bundled Simplification of Trail Sets[END_REF] for bundled simplification of 2D and 3D curve-sets. New curves are next generated after deforming the cluster centroids to manipulate the underlying data in a simple way while preserving their statistical properties. The problem of short to mid-term aircraft trajectory prediction is considered in [START_REF] Tastambekov | Aircraft trajectory forecasting using local functional regression in Sobolev space[END_REF]. This approach is based on a local linear functional regression using wavelet decomposition. In [START_REF] Gregorutti | Grouped variable importance with random forests and application to multiple functional data analysis[END_REF], random forests for functional data are used for minimizing the risk of accidents and identifying explanatory factors in the context of aviation safety. In the closed context of road traffic, the derivative information was used in a functional analysis of speed profiles of vehicles [START_REF] Andrieu | Functional modeling of speed profiles adapted to the infrastructure and methodology of construction of an aggregated speed profile[END_REF], providing a functional modeling of space-speed profiles of vehicles that leads to a complex nonparametric regression problem.

A major asset of working with functional data framework is the ease of adding a priori information by carefully selecting the Hilbert space and basis. 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.

Moreover, the choice of the L p metric as a measure of similarity between trajectories may be problematic in many situations and does not reflect the high internal structure of aircraft trajectories. As noted in [START_REF] Srivastava | Functional and Shape Data Analysis[END_REF], this metric suffers from several drawbacks, especially in registration of functional data when functions are shifted, owing to time lags or general differences in dynamics and alternative metrics should be implemented in a functional data analysis.

For overcoming these difficulties, several methods are based on warping distance such as DTW [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF], LCSS [START_REF] Vlachos | Discovering similar multidimensional trajectories[END_REF], EDR [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF] and ERP [START_REF] Chen | On the marriage of lp-norms and edit distance[END_REF]. Alternative distances focus on the shape of trajectories, such as the Hausdorff distance [START_REF] Hausdorff | Grundz uge der mengenlehre[END_REF], the Fréchet distance [START_REF] Fréchet | Sur quelques points du calcul fonctionnel[END_REF] and the SSPD distance [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF]. A comprehensive review of these distances used in distance-based trajectory clustering is given in [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF].

B. The Riemannian framework

An alternative solution consists in switching from a functional data framework to a geometrical one, by representing trajectories in a Riemannian framework. Indeed, curves can be made amenable to functional data statistics by representing them as points on the so-called shape manifold, that is formally defined as a quotient of the manifold of smooth immersions from [0, 1] to R n . It is not a Hilbert space, but it can be provided with Riemannian metrics that allow geodesic distance computation, although some care must be taken to avoid degeneracy of the metric [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF].

In the shape space literature, curves are considered as geometric objects and are represented in the so-called shape space, that is the set of immersions quotiented by the group of smooth diffeomorphisms of the interval [0, 1]. An overview of various notions of shape spaces is given in [START_REF] Bauer | Overview of the Geometries of Shape Spaces and Diffeomorphism Groups[END_REF]. Three family of metrics on shape space can be distinguished: the almost local metrics, the Sobolev-type metrics and the right-invariant metrics on the diffeomorphism group of the ambient space. Note that this kind of metrics verify the reparametrization invariance property.

More general almost local metrics and Sobolev-type inner metrics on the space of plane curves were considered in [START_REF] Michor | An overview of the riemannian metrics on spaces of curves using the hamiltonian approach[END_REF] and some algorithms [START_REF] Mumford | The geometry and curvature of shape spaces[END_REF] use a shape manifold representation in order to derive a metric between sample paths. The special case of order one Sobolev type metric is presented in [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF] and a similar approach takes into account a greater amount of information on the space separating two curves by computing the distance in the manifold itself rather than in one tangent plane [START_REF] Brigant | Computing distances and geodesics between manifoldvalued curves in the SRV framework[END_REF].

Having it at hand, statistical analysis of functional data such as trajectories clustering may be performed in a standard way. It worth mention that most of the time a mean of computing the centroid of a set of curves is mandatory. In [START_REF] Srivastava | Functional and Shape Data Analysis[END_REF], elastic metric for curves has been integrated into statistical analysis of functional data. This so-called Elastic Functional Data Analysis performs alignment, principal component analysis, and modeling of multidimensional and unidimensional functions using the square-root velocity framework developed in [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF], [START_REF] Tucker | Generative models for functional data using phase and amplitude separation[END_REF]. Nice libraries such as the R package fdasrvf have been developed [START_REF]CRAN package fdasrvf[END_REF].

Considering methods falling in this category, it appears that the requirements of the landing trajectory analysis are not fulfilled: enforcing a full parametrization invariance as in the case of shape-based approaches prevents the use of longitudinal acceleration that enters the non-slip condition. On the other side, Sobolev-like metrics are designed to be tailored to specific needs. Motivated by a real use case where one wants to assess runway adherence condition by observing only the radar tracks of landing aircraft, we propose in this paper to use a Riemannian framework for curves with velocity information.

C. Related works in Trajectory Clustering

The mathematical problem arising from the estimation of the runway adherence condition falls within the general framework of curve clustering. Several well established algorithms may be used for performing clustering on a set of trajectories.

First, in sample-based approaches, time samples on trajectories are used in a classical multivariate analysis. Due to the huge number of points, a dimension reduction procedure must be applied first. Random projections followed by a spectral clustering was successfully applied in a study conducted by the Mitre corporation on behalf of the Federal Aviation Authority (FAA) [START_REF] Enriquez | Identifying temporally persistent flows in the terminal airspace via spectral clustering[END_REF]. The most important limitation of this approach is that the shape of the trajectories is not taken into account when applying the clustering procedure unless a resampling procedure based on arclength is applied.

Secondly, density-based algorithms rely on a density map on the state space and a definition of similarity measure. They aim to attract trajectories towards areas of high density and iterate the process until convergence. This kind of algorithms are well suited for already structured graph structure, most implementation derived from the original work presented in [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. The distance-based density-based algorithms developed in [START_REF] Mahrsi | Graph-based approaches to clustering network-constrained trajectory data[END_REF], [START_REF] Kim | Spatial and temporal characterization of travel patterns in a traffic network using vehicle trajectories[END_REF] are well adapted to road traffic as vehicles are bound to follow predetermined segments but are unable to work on an unknown airspace organization.

Due to the functional nature of trajectories, that are basically mappings defined on a time interval, it seems more appropriate to resort to techniques based on times series as surveyed in [START_REF] Liao | Clustering of time series data -a survey[END_REF], [START_REF] Rani | Recent techniques of clustering of time series data: A survey[END_REF], or functional data statistics, with standard references [START_REF] Ramsay | Functional Data Analysis[END_REF], [START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice, ser. Springer Series in Statistics[END_REF]. Again, in both approaches, a distance between pairs of trajectories or, in a weaker form, a measure of similarity must be available. The algorithms of the first category are based on sequences, possibly in conjunction with dynamic time warping [START_REF] Meesrikamolkul | Shape-based clustering for time series data[END_REF], while in functional data analysis, samples are assumed to come from an unknown underlying function belonging to a given Hilbert space. Both approaches yield similar end algorithms, since functional data revert for implementation to usual finite dimensional vectors of expansion coefficients on a suitable truncated basis. For the same reason, model based clustering may be used in the context of functional data even if no notion of probability density exists in the original infinite dimensional Hilbert space as mentioned in [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. A nice example of a model-based approach working on functional data is the R package funHDDC [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF], [START_REF]CRAN package funHDDC[END_REF].

The purpose of this paper is to introduce a new kind of Riemannian metric that is especially adapted to the clustering of curves where the velocity is a discriminating feature. In such a case, the original shape space approach cannot be used since the parametrization invariance will wipe out the velocity information. A partial parametrization invariance is introduced, yielding a bundle shape space model on which a relevant metric can be defined. The design of the metric was based on the equations of motion and reflects the internal structure of the data. The numerical implementation of the path computation of the new metric is given. For simplicity, we will use a discrete approximation of the total energy of the sampled path rather than a practical computation algorithm based on a shooting method. The minimal approximate energy will be obtained by a standard optimization algorithm and will be next used in the clustering phase.

The paper is organized as follows. In Section II, we formally describe the theoretical framework and define an adapted metric taking into account the velocity information. In Section III, the numerical implementation of the new metric and the K-medoids clustering method [START_REF] Kaufman | Clustering by Means of Medoids[END_REF], [START_REF]Finding Groups in Data: an introduction to cluster analysis[END_REF] is described in a general setting. This algorithm allows to implement another metric than the usual Euclidean distance. Next, Section IV illustrates the performance of clustering algorithms with the new metric for low adherence detection on a set of simulated and real aircraft trajectories. These results are also compared with competing distances. Finally, some comments and future works are drawn in Section V.

II. THEORETICAL FRAMEWORK

This section is dedicated to the definition of a relevant metric for curves taking into account velocity information. The framework presented here will overcome the limitations of the shape approach through the use of a metric in which the slip characteristic of the landing path is explicitly considered.

A. Slip detection

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 high braking action or a 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 ) (1)
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 (2)
with M the aircraft mass, g the gravity of 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] Ramsay | Functional Data Analysis[END_REF]. At that point, slip occurs and µ drops to the value of the dynamic friction coefficient µ d . The contact force 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 expression µ as a function of λ [START_REF] Rajamani | Vehicle Dynamics and Control, ser. Mechanical Engineering Series[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 ( 3
)
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 (3) 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 (4)
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 can be deduced from the previous remark by considering the angle θ between the acceleration and speed vectors:

sin(θ(s)) = κ(s) D s γ(s) 2 D ss γ(s) = det (D s γ(s), D ss γ(s)) D s γ(s) D ss γ(s) (5)
where γ is the aircraft trajectory and κ its curvature. In the above expression and through all the document, the symbol D s stands for the partial derivative with respect to variable s. Higher order derivatives are written similarly as D s1...s1,s2...s2,... by repeating the variable p times to indicate a partial derivative of order p.

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 s γ, D ss γ). Reciprocally, under slip conditions, a trade off has to be made between path following and deceleration: the angle θ will thus increases towards the limiting value ±π/2.

Please note that θ is not defined in the portions of the trajectory where the acceleration or the velocity vanishes. While the second situation is highly uncommon, as it will indicate a stop during landing or taxiing, the first one can be encountered when the aircraft is not braking nor turning. In such parts of the trajectory, it is not possible to infer an adherence condition, and curve comparison can be done only on a geometrical basis. In the addressed application, aircraft in the observed part of the landing trajectory are decelerating, so that the above issue is not a concern.

B. An adapted metric in the space of trajectories

Recall that a smooth 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 ). It has the structure of Banach manifold, its tangent space at γ ∈ Imm([0, 1], R 2 ) being the vector space C ∞ ([0, 1], R 2 ). A tangent vector is thus a couple (γ, h) with γ the base curve and h an element of C ∞ ([0, 1], R 2 ) that is interpreted as an infinitesimal displacement field along γ. This point can be formalized using the following notion of admissible variation. Definition 1. Let γ be a smooth loop. An admissible variation of γ is a smooth mapping

Φ : ]-, [×[0, 1] → R 2 , > 0, such that Φ(0, •) = γ(•).
An admissible variation defines a tangent vector (γ, D t Φ| t=0 ). 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 is used to derive a Riemannian metric on the quotient space Imm(S 1 , R 2 )/Diff ([0, 1], R 2 ). These notions can be easily generalized to the case of open curves.

In the present work, a similar approach will be taken. However, due to the fact that the slip condition must come into play, it is not meaningful to keep invariance under change of parametrization. Instead, a weaker invariance by affine change of parametrization will be obtained at the end. Furthermore, curves with vanishing second derivative must be excluded since the slip angle θ in [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF] is not defined at points where D ss γ(s) = 0. The last condition boils down to the requirement that the curve s ∈ [0, 1] → (γ, D s γ) be an immersion. The space of such objects will be denoted by Imm * ([0, 1], R 2 ). An mentioned above, this is not an issue for the application, due to the restriction of the braking part of the trajectory. In a more general setting, it may be necessary to segment curves into non vanishing second derivative parts and compare only them.

The slip angle in (5) has a nice variational interpretation as indicated in the next lemma.

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

D t A(0) = 1 0 det (D t Φ(0, s), D s Φ(0, s)) ds (6)
where A(t) is the net area enclosed by the loop

Φ(t, •), φ for t ∈] -, [.
Proof. For any t ∈] -, [, the loop Φ(t, •), φ is piecewise smooth so that the net area A(t) can be computed using Stoke's formula:

A(t) = 1 2 1 0 det(Φ(t, s), D s Φ(t, s))ds (7)
+ 1 2 1 0 det(φ(s), D s φ(s))ds (8)
Taking the derivative with respect to t yields:

D t A(t) = 1 2 1 0 det(D t Φ(t, s), D s Φ(t, s))ds (9)
+ 1 0 det(Φ(t, s), D st Φ(t, s))ds (10)
Recalling that the end points of Φ(t, •) are fixed, one can use integration by parts to get:

D t A(t) = 1 2 1 0 det(D t Φ(t, s), D s Φ(t, s))ds (11)
-

1 0 det(D s Φ(t, s), D t Φ(t, s))ds (12)
and finally, using the alternating property of the determinant:

D t A(t) = 1 0 det(D t Φ(t, s), D s Φ(t, s))ds (13)
and the conclusion follows by letting t = 0.

Using Lemma 1, the integral :

1 0 |sin(θ(s))| D s γ(s) ds = (14) 1 0 |det (D s γ(s), D ss γ(s))| D ss γ(s) ds (15)
may be interpreted as the total infinitesimal area swept by the curve γ when moved in the direction D ss γ.

The integral in ( 14) is invariant under affine change of parametrization. Given a curve γ defined on an arbitrary interval [a, b] , it is thus possible to go back to fixed interval [0, 1]. As a consequence, any curves will be assumed to be defined on [0, 1]. Please note however that general invariance under change of parametrization cannot be obtained, and will not be meaningful as the information sought after is clearly dependent on the velocity.

The slip angle θ can be obtained in a convenient way using the fact that the curves take their values in R 2 and can thus also be considered as complex valued mappings. Within this frame, it comes for the expression of the slip angle θ if curve γ:

e iθ = D s γ |D s γ| D ss γ |D ss γ| . ( 16
)
Using the complex representation, the following proposition gives access to the variation of the slip angle.

Proposition 1. Let Φ : ] -, [×[0, 1]
to C be an admissible variation of γ. Let θ(t, s) be the slip angle at s of the curve s → Φ(t, s). Then:

D t e iθ = 1 |D s Φ| D st Φ| N D ss Φ |D ss Φ| N + 1 |D ss Φ| D sst Φ| N D s Φ |D s Φ| N
where the notations N and | N (resp. N and | N ) denotes the unit normal vector and the component normal to the curve

s → Φ(t, s) (resp. s → D s Φ(t, s)).
Proof. The proposition is a direct computation of the derivative. The only point to note is that:

D t D s Φ |D s Φ| = D st Φ |D s Φ| - D s Φ |D s Φ| D st Φ.D s Φ |D s Φ| 2 (17) = 1 |D s Φ| D st Φ| N N (18)
and the equivalent derivation for D ss Φ:

D t D ss Φ |D ss Φ| = 1 |D ss Φ| D sts Φ| N N .
Proposition 2. Under the same assumptions as in Proposition 1:

(D t θ) 2 = D st Φ| N |D s Φ| - D sst Φ| N |D ss Φ| 2
Proof. Using the fact that |e iθ | = 1 and Proposition 1, it comes:

(D t θ) 2 = 1 |D s Φ| 2 |D st Φ| N | 2 + 1 |D ss Φ| 2 |D sst Φ| N | 2 + 2 1 |D s Φ||D ss Φ| D st Φ| N . T .D sst Φ| N T with: T = D s Φ |D s Φ| , T = D ss Φ |D ss Φ| .
Since the curve is planar, N = iT (resp. N = i T ). This combined with the expression of the normal component

D st Φ| N = D st Φ, N (resp. D sst Φ| N = D sst Φ| N , N ) yields: D st Φ| N . T .D ss D t Φ| Ñ T = -D st Φ, N D sst Φ| N , N
and the result follows.

An admissible variation Φ of a curve γ defines a tangent vector to γ in the space Imm * ([0, 1], R 2 ) by D t Φ| t=0 . Such a tangent vector is thus a couple (γ, u) where u ∈ C ∞ ([0, 1], R 2 ). The above expression for (D t θ) 2 allows to define a semi-norm in the tangent space to the curve γ by integration along the curve. Definition 2. Let (γ, u) be a tangent vector to γ. The seminorm M (γ,u) is defined by:

M 2 (γ,u) = 1 0 D s u| N |D s γ| - D ss u| N |D ss γ| 2 ds
with N (resp. N ) the normal vector to γ (resp. D s γ).

Given two curves γ 1 , γ 2 , a path between them in

Imm * ([0, 1], R 2 ) is smooth mapping Φ : [0, 1] × [0, 1] → R 2 such that: • Φ(0, •) = γ 1 (•), Φ(1, •) = γ 2 (•) • For all t in ]0, 1[, Φ(t, •) ∈ Imm * ([0, 1], R 2 ).
Definition 3. Let Φ be a path between γ 1 and γ 2 . Its θdeformation energy is defined as:

E(Φ) = 1 0 |D t Φ(t, 0)| 2 + |D t Φ(t, 1)| 2 dt + 1 0 1 0 D st Φ| N |D s Φ| - D sst Φ| N |D ss Φ| 2 ds dt.
The first term is needed as the second part involves only derivatives of the curve and will not take endpoints into account. The above definition of the energy is valid for curves in Imm * ([0, 1], R 2 ), but is not invariant by affine change of parametrization. Furthermore, it is based on a semi-metric, resulting in possible degeneracies. The path energy is thus modified according to Definition 4.

Definition 4. Let Φ be a path between γ 1 and γ 2 . Its total deformation energy is defined as:

E(Φ) = 1 0 |D t Φ(t, 0)| 2 + |D t Φ(t, 1)| 2 dt + 1 0 1 0 |D t Φ| 2 |D s Φ| ds dt + λ 1 0 1 0 D st Φ| N |D s Φ| - D ss D t Φ| N |D ss Φ| 2 |D s Φ| ds dt
where λ > 0 tunes the relative importance of shape and slip angle.

The first double integral is a Sobolev Riemannian metric that accounts for shape deformation. The free parameter λ must be tuned in applications. The energy allows the computation of geodesic path between curves.

Definition 5. The geodesic path Φ 0 between γ 1 and γ 2 is, when it exists, defined by:

Φ 0 = argmin Φ∈I E(Φ)
where I is the set of paths from γ 1 to γ 2 .

The energy of a geodesic path defines a similarity measure between curves. It will be used in the sequel within a clustering algorithm.

III. NUMERICAL IMPLEMENTATION

This section describes the numerical implementation of the path computation of the new metric described in Section II-B. Then, the clustering procedures will be presented.

A. Geodesic path computation

Two approaches can be used for computing a geodesic path Φ 0 minimizing the total energy. The first one relies on a system of differential equations involving D t Φ that is solved starting at t = 0 from the first curve. An iterative shooting procedure reduces gradually the distance at t = 1 to the second curve, yielding an approximate geodesic path. This algorithm is appealing as it does not require an approximation to the entire path Φ 0 , but just of the curve s → Φ(t, s) for a fixed t. Its two main drawbacks are first the need to establish the equation of geodesics, that may be quite complicated, and second the numerical instabilities in the shooting phase. It will nevertheless be investigated in a future work. The second approach is to approximate the full path Φ by its samples Φ ij = Φ(x ij ), i = 0 . . . n, j = 0 . . . m on a grid defined by the points x ij = (t j , s i ). In the current implementation, the grid is evenly spaced. A numerical differentiation procedure is applied to obtain the derivatives D s , D ss , D t and takes the form of matrices, still denoted with D s , D ss , D t , that are applied to the matrix M = (Φ ij ) of path sample points. Since the derivative with respect to s must be applied columnwise, the derivative D s Φ is approximated by D s .M . Conversely, derivation is performed rowwise for D t Φ, so that the approximation is M.D t . Finally, the first and the last column of Φ are held constant, since they represent sampled version of the endpoint curves γ 1 , γ 2 .

Using the discrete approximation, the total energy of the sampled path described by the matrix M is a sum of three terms, hereafter denoted by E 1 , E 2 , E 3 . Their respective expressions are given in equations

E 1 (M ) = m j=0 |(M.D t ) 0j | 2 + |(M.D t ) nj | 2 (19)
E 2 (M ) = m j=0 n i=0 |(M.D t ) ij | 2 |(D s .M ) ij | (20)
E 3 (M ) =λ m j=0 n i=0 (D s .M.D t | N ) ij |(D s .M ) ij | (21)
-

(D ss .M.D t | N ) ij |(D ss .M ) ij | 2 |(D s .M ) ij |. ( 22
)
The optimal approximate path is obtained by a standard numerical optimization algorithm, chosen here to be LM-BFGS [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF] that is well suited to high dimensional problem. Initializing Φ with a linear homotopy ensure convergence in a few iterations on most problems. The minimal energy obtained is used in the following clustering phase.

B. Clustering methods

K-medoids algorithm [START_REF] Kaufman | Clustering by Means of Medoids[END_REF], [START_REF]Finding Groups in Data: an introduction to cluster analysis[END_REF] is an efficient clustering technique that aims to partition unlabelled data into clusters. Unlike the well known K-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] which is sensitive to outliers, K-medoids algorithm is more robust and uses an actual point in the cluster to represent it. As this method produces realisctic cluster centers, it may be preferred for aircraft trajectory clustering. The most benefit of K-medoids algorithm compared to the K-means is that there is no need to explicitly compute the cluster centers since we just need to store the indices of the trajectories minimizing the total distance to other trajectories in that cluster. Note that K-medoids requires the use of a metric distance that satisfies the triangular inequality. It may be a serious drawback for distances that are not metrics. As recommanded in [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF], a hierarchical cluster analysis (HCA) should be prefered.

As a matter of fact, it is possible to perform the clustering with a pairwise distance matrix as input data instead of the raw data. Thereby, we can derive an algorithm with the metric defined in Section II in order to classify a set of n spatiotemporal trajectories {x 1 , . . . , x n }. Similarly to K-means and other clustering techniques such as HCA, the number of clusters K is assumed to be known. In practice, the knowledge of the data can help the user to estimate the value of K. In the case of aviation safety, one could choose K = 3 corresponding to three levels of alerts: low, medium, high. A fourth cluster may be added in the procedure for detecting outliers trajectories. In general, the number of clusters is not known and has to be computed from the data. In [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF], the authors introduced the Silhouette criterion as a graphical method for interpreting and validating the results of a clustering algorithm. Broadly speaking, the Silhouette score measures the tightness and the separation of the clusters. More formally, for an observed trajectory x i , it is defined as:

s(i) = b(i) -a(i) max(a(i), b(i))
,

where a(i) is the distance between x i and all other trajectories in the same cluster and b(i) is the distance between x i and all other trajectories in the next nearest cluster. The overall Silhouette score, which ranges between -1 and 1, is then the mean value of the score computed on each sample, i.e.

s = 1 n n i=1 s(i).
Hence, the number of clusters K can be chosen by maximizing the Silhouette criterion. It can also be used as a validation metric of a clustering method: the higher value of s is, the better is the clustering. Note that the mean values of a(i) and b(i) can be interpreted respectively as the intra-cluster variance and the inter-cluster variance. Knowing that, a higher value of s indicates that the clusters are dense and well separated while a lower value of s refers to an incorrect clustering. A score around zero indicates overlapping clusters. Note also that the K-medoids can be sensitive to the initialization of the algorithm and may find a local minimum. This is actually a classical drawback of the clustering methods. A common solution is to repeat the clustering with different initializations and to retain the best one with respect to the Silhouette score.

IV. NUMERICAL EXPERIMENTS In the following, the new metric defined in Section II is called the slippery metric. The slippery metric described is tailored to take into account information on acceleration and should be able to detect abnormal deceleration levels that may indicate a potential degraded runway condition. This metric has an affine invariance property by change of parametrization, speed and acceleration of trajectories remaining unchanged except for a multiplicative constant. In this section, the performance of this metric is compared with other existing competitors: one geometric metric, the Square Root Velocity (SRV) metric [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF]- [START_REF] Tucker | Generative models for functional data using phase and amplitude separation[END_REF], and seven other distances described in [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF].

Geometric metrics are relevant for clustering trajectories having different morphologies, and for which second-order information is not required; this is the case when some aircraft trajectories deviate from the nominal one. As noted in [START_REF] Brigant | Geomstats: Shape analysis of curves with the Square Root Velocity metric[END_REF], the SRV metric is invariant under change of parametrization: if the two curves are reparametrized in the same way, the distance does not change. However, if the curves are reparametrized in different ways, the distance changes. Even if the SRV metric is not tailored to take into account information on acceleration, it might be able to detect parametrization differences between curves in addition to changes in the shape. Usually, an optimal matching procedure is used to produce a distance, that can detect differences in the shapes of the curves irrespective of their parametrization. This should not take into account the velocity information and is not be suitable for slip detection. In this section, we will consider both kinds of SRV metric: the first one will be called the SRV metric before alignment and the second one will be called the SRV metric after alignment. We have used the algorithm developed in [START_REF] Brigant | A discrete framework to find the optimal matching between manifold-valued curves[END_REF], implemented in the Geomstats python package [START_REF] Brigant | Geomstats: Shape analysis of curves with the Square Root Velocity metric[END_REF], [START_REF]Geomstats: open-source Python package Geomstat for computations, statistics, and machine learning on nonlinear manifolds[END_REF].

In addition to the SRV metrics, [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF] provides a comprehensive review of different distances used in the literature, especially in the context of vehicle trajectory clustering. But most of the used competitor distances are not metrics. For this reason, in the following first subsection, we will perform a hierarchical cluster analysis (HCA) in a validating step on simulated datasets. In the second subsection, the performance of the two SRV metrics and the slippery metric will be compared by using a K-medoids clustering procedure to a real dataset of landing aircraft tracks.

A. The simulation procedure

As already mentioned, it is virtually impossible to ensure that an observed trajectory that deviates from the nominal one is the result of a bad adherence condition. To perform such a task, a physical measurement is mandatory and the only way to correlate runway quality to landing tracks is to have recorded trajectories right after or before a measurement.

A more tractable approach is to simulate landing and taxiing in a realistic fashion, so that the detection performance of the methods can be estimated. Very good models of wheeled ground vehicles are available in the literature [START_REF] Rajamani | Vehicle Dynamics and Control, ser. Mechanical Engineering Series[END_REF] and the case of braking aircraft has been considered in [START_REF] Jones | Modélisation des forces de contact entre le pneu d'un avion et la piste[END_REF]. For the purpose of the study, a very detailed model is not needed, but a good conformance to observations is required. The approach taken was to use the model of Figure 1. The aircraft is controlled Fig. 1: Aircraft physical model. using the engine thrust and the front wheel steering angle.

landing gear tires train is replaced by a single equivalent tire. The simulation is performed by solving the equations of motion using a Runge-Kutta 45 scheme with adaptive step. The control loop consists of two Proportional-integratedderivative (PID) regulators: the inner acts on the steering angle of the front wheel and corrects lateral deviations from the intended path while the outer controls deceleration. When the mobile is skating, the first control loop takes precedence and consumes the apart of the available friction force to minimize the lateral error. The remaining friction is transfered to the outer loop to apply some deceleration. It is clear that this process will increase the slip angle in a trajectory turn.

A simulator was coded in Java and applied to synthetic tracks with tunable adherence coefficient. Since slippery can be detected easily within the simulation loop, all output trajectories are flagged as "skating" or "normal". Two datasets were generated from this simulator to compare the slippery metric with all others competitors. In order to more precisely compare the behavior of the SRV metrics and the slippery metrics, other more tractable datasets were simulated following a circle shape and different parametrizations of curves.

B. The aircraft trajectory simulator datasets and comparison with all competitors 1) Simulated trajectories: Landing and taxiing trajectories are simulated following a circle representing the nominal trajectory. The initial landing point is randomly disturbed as well as the initial velocity. In this way, it may induce a skating at the beginning of the landing trajectory which is then corrected by the regulators as described above. Note that, following the action of the PID regulators, simulated skating trajectories end before the normal ones. Moreover, all the paths take the same initial direction following the nominal trajectory as illustrated in Figure 2. The dataset of simulated trajectories is composed of 16 skating trajectories and 284 normal trajectories. As the analysis of landing trajectories requires to take into account information on acceleration, simulated trajectories have not been normalized. The slippery metric is compared with both SRV metrics and the distances described in [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF]: some warping-based distances such as DTW [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF], LCSS [START_REF] Vlachos | Discovering similar multidimensional trajectories[END_REF], EDR [START_REF] Chen | Robust and fast similarity search for moving object trajectories[END_REF] and ERP [START_REF] Chen | On the marriage of lp-norms and edit distance[END_REF], and some shape-based distances, such as the SSPD distance [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF], the Hausdorff distance [START_REF] Hausdorff | Grundz uge der mengenlehre[END_REF] and the Fréchet distance [START_REF] Fréchet | Sur quelques points du calcul fonctionnel[END_REF]. The discrete Fréchet distance was preferred, as requiring a lower computational cost. Because the discrete Fréchet distance, just like SSPD, DTW, LCSS and EDR distances are not metric, a HCA clustering method is then performed in order to detect two clusters: the normal trajectories cluster and the skating trajectories one. In this framework, complete, median, ward and centroid linkage criteria have been compared. Finally, following [START_REF] Besse | Review and perspective for distance-based clustering of vehicle trajectories[END_REF], the ward criterium has been chosen, usually giving the best results and producing spherical clusters. Note that the LCSS and EDR warping-based distances require to set a spatial threshold that makes the distance robust to noise. A high threshold value can produce all similar distances, while by using a too low threshold, only very close curves are considered similar. We set this threshold to 0.01 after comparing different results. Table I displays the number of trajectories in each cluster returned by the HCA procedure. The output trajectories are then flagged as "skating" or "normal". The slippery metric and both SRV metrics, as well the ERP warping-based metric, the Hausdorff and discrete Fréchet shape-based distances, perform well and correctly detect the 16 skating trajectories cluster and the 284 normal trajectories one. The DTW warping-based distance produces one cluster containing only normal trajectories and one cluster composed of all the skating trajectories, but also a large number of normal trajectories. The EDR and LCSS warping-based distances and the SSPD distance detect one cluster with some skating trajectories, the others being clustered with the normal trajectories. Table II displays some evaluation measures for DTW, EDR, LCSS and SSPD distances: the false negative rate (FNR), the false positive rate (FPR) and some indices measuring how similar the returned clusters are to the skating and normal labels. FNR represents the miss rate, that is the rate of skating trajectories in the cluster of returned normal trajectories corresponding to an undetected skidding, distances with high value miss rate are not suitable from a safety point of view. FPR corresponds to the false alarm rate, that is the rate of normal trajectories falsely detected as skating. A high FDR value corresponds to a situation generating significant costs for ATM as noted in Section I. The Matthews correlation coefficient takes on a value between 0 and 1, with 0 meaning that the two returned clusters are completely different from the skating and normal clusters. Higher values indicate greater similarity and an indicator of 1 means that the output trajectories are correctly clustered as skating and normal. Even if the DTW miss rate is 0, the DTW false alarm rate is very high (92%), as the SSPD distance (90%). EDR and LCSS both produce the lowest miss rate and false alarm rate values, but the Matthews correlation coefficient indicates that more than half of skidding trajectories are not detected by using EDR and LCSS distances.

As noted previously, skating trajectories can be distinguished from normal ones due to their end points. This shape argument could explain why the SRV metrics and the Hausdorff, discrete Fréchet and ERP distances work well for clustering skid and non-skid situations. For this reason, we have truncated the simulated trajectories, so that all trajectories end almost at the same location. We assess clustering for ERP, discrete Fréchet, Hausdorff distances, both SRV metrics and the slippery one in Table III Results in Table III show that the ERP metric, as well as the Hausdorff and Discrete Fréchet distances, no longer correctly detect the two skating and normal clusters. the low miss rate value shows that few skating trajectories are falsely detected as normal, the false alarm rate value is very high. Moreover, the Matthews correlation coefficient indicates that the predicted clusters are not similar to the skating and normal clusters.

Finally, only the slippery metric and both SRV metrics seem to be able to correctly cluster trajectories in a skid and nonskid situation. This result is surprising, particularly that of the SRV with alignment, which is probably explained by a weak numerical bias, detected by the HCA procedure. The result from the SRV without alignment is coherent with the fact that, before the optimal matching procedure, the differences in parametrization may be detected, when differences in shapes are especially not enough significant to hide this effect.

2) Slightly deformed simulated trajectories: In order to assess the robustness of the results obtained by using the slippery metric and the SRV metrics, we have slightly deformed trajectories by using a random coefficient (around 1%) so that simulated landing trajectories are now randomly flattened following an ellipsoid and are more realistic.

When HCA is performed on the entire trajectories, both SRV metrics and the slippery metric perform well. We then focused on the beginning of the trajectories before the PID regulators act skidding. Now, both SRV metrics can no longer correctly detect the two clusters while the slippery metric can do so. Results from Table V show that almost all skating trajectories are clustered with normal ones. These two metrics produce a very high false alarm rate (around 90%) and the low Matthews correlation coefficient indicates that the clustering algorithm does not correctly distribute the trajectories into clusters. In Figure 3, we can visualize the clusters produced

Fig. 3: HCA with (a) the slippery metric, (b) the SRV metric without alignment and (c) the SRV metric with alignment. Two predicted clusters: cluster 1 (blue) and cluster 2 (red).

from HCA with the slippery metric and both SRV metrics. When using the SRV metrics, almost all the skating trajectories are clustered with the normal trajectories whose the starting points have low values. It seems that normal trajectories are distributed in the two predicted clusters in equal proportions relatively to the starting point, the skating trajectories being clustered in one of these clusters.

Finally, it seems that, among all the competitors, only the SRV metrics might be able to correctly detect the skidding in certain situations. These results are quite surprising in the case of deformed simulated trajectories. We suspect that the differences in the shapes along the entire trajectories are not sufficiently high to hide the parametrization differences or a numerical bias. However, the datasets generated by the simulator may eventually suffer from biases due to the action of the PID regulators. This is the case of the end points of trajectories, but other effects may be present.

C. Application to more tractable simulated datasets

In order to highlight the effect of shape differences, we have compared results from the SRV metrics and the slippery metric with more tractable datasets. In the following datasets, we have considered two profiles of velocity along the curves: a constant speed profile corresponding to skating and a braking one at the beginning of the curves. The datasets of simulated curves are composed of 30 curves, mixing the two profiles of velocity in equal proportions.

1) Trajectories of the same shape: We first consider trajectories following the shape of a circle with initial point being randomly disturbed. Each curve i is a set of positions

traj(i) = (a i cos θ + c i , b i sin θ + c i ), θ ∈ [0, π/2],
where c i is evenly distributed in an interval of amplitude α and a i = b i = 1. The SRV metric is not sensitive to differences in tranlation and should vanish. Surprisingly, the HCA procedure can detect the two clusters of speed profile, even in the case of the SRV without alignment. When examining the distance matrices produced by the SRV metric with alignment, the distance values of similar speed profile curves is around of 10e-14 and 10e-5 for different speed profile curves. These differences seem negligible at first glance, however they are anyway detected by a HCA procedure. In this case, the clustering result is clearly due to a numerical bias. On the other hand, the SRV metric without alignment produces similar distance values for similar speed profile curves, but non negligible values for different speed profile curves around 0.15. This result is less surprising because, as noted in [START_REF] Brigant | Geomstats: Shape analysis of curves with the Square Root Velocity metric[END_REF], the SRV metric may be sensitive to differences in reparametrization.

2) Trajectories with deviations from the initial shape: We consider the circle-shaped trajectories as in the previous example. However, some of them exhibit a Gaussian-like deviation from the nominal circle-shaped curve. This corresponds to the shape of landing trajectories with lateral deviations from the centreline. For each curve i, this deviation d i can be inflated by a multiplicative coefficient r as follows

d i (t) = (r * exp(-0.5 * (t -0.5) 2 /σ 2 ), t ∈ [0, 1],
where t represents a time location and σ a fixed standard deviation. Two situations were considered with different values of the coefficient r : low-medium deviations (r = 1%, 5% and 10%) and medium-high deviations (r = 1%, 10% and 30%). Again, for each scenario, the two speed profiles are equally distributed into the three subgroups curves. Table VII shows that, with low-medium deviations, the SRVF metric with alignment can no longer correctly detect the two clusters. Both skating and braking trajectories are similarly clustered into the predicted clusters in equal proportions. Moreover, the miss rate (50%), the false alarm rate (50%) and the zero Matthews correlation coefficient mean that the clustering algorithm randomly distributes the trajectories into clusters. Note that the SRVF metric without alignment can detect skating and braking when deviations are relatively low. Table VIII displays the results when the differences in the shape are higher. The SRV metric without alignment now favors shape detection to the detriment of speed profile detection. In Figure 6, we can visualize that both SRV metrics cluster the curves relatively to the shape differences while the slippery metric performs well. 3) Trajectories with slightly deformed deviations: Finally, the previous datasets are slightly deformed by adding Gaussian distributed noise N (0, ε = 0.05) along the curves. A smoothing splines regression technique is next applied to recover smooth trajectories. The resulting dataset is then composed of slightly deformed curves relatively to the previous example. Both SRV metrics are now unable to correctly cluster the curves as skating and braking. Even if the evaluation criteria are better for the SRV without alignment than for the SRV with alignment, the false alarm rate remains important (20%). Finally, with medium-high deviations, we obtain the same results as those found in the previous examples: the parametrization detection is hidden by the shape differences by both SRV metrics while it is clearly identified by the slippery metric. These results highlight those obtained with the aircraft trajectory simulator datasets. When the curves have the same shape, both SRV metrics perform well, due to a numerical bias for the SRV metric with alignment. When the curves have different morphologies, the SRV with alignment behaves as expected and detects only differences in shape. The numerical bias seems to be hidden in this case. This observation does not hold for the SRV without alignment, which depends on the parametrization of the curves. Even if this metric is not tailored to take into account information on acceleration, it is able to perform well in certain situations when shape differences are not high enough to hide differences in parametrization. As an alternative, the slippery metric that is tailored to detect the information on acceleration performs well in each studied case.

D. Application to Radar tracks

In this section, we apply the K-medoids clustering method to a database of radar tracks in order to detect bad runway condition and skidding. Note that an HCA procedure was also performed and analogous results were obtained.

1) Description of the data: The Radar tracks used here are a set of landings coming from A-SMGCS 1 which is a system that provides air traffic controllers with services to maintain airport safety and capacity. These data are generated by sensor fusion of Primary Surveillance Radar (PSR), Secondary Surveillance Radar (SSR) using onboard transponders and multilateration Automatic Dependent Surveillance Broadcast (ADS-B).

Each landing trajectory T is defined as a set of spatiotemporal data, i.e. includes GPS position x k , y k and the aircraft speed v k :

T = (p 1 , . . . , p n ) = (γ(t 1 ), . . . , γ(t n )),
where p k = (x k , y k , v k ) and γ : [0, 1] -→ R 3 . Note that deceleration values are derived from aircraft speed as a preprocessing step. The dataset consists of 357 landing trajectories.

We have chosen to focus on the runway exit curve only since slippery is mainly visible in curves. Then, following experts, as skidding is a loss of grip with the runway, a degraded runway condition is defined as follows:

• The aircraft deviates from the centreline • The deceleration value is greater than the runway grip coefficient (which is low when the runway is slippery) Following this definition, skidding is not only reduced to a lateral deviation of the trajectory. The aircraft deceleration value is indeed used as the indicator of the runway state. A low deceleration level indicates that the runway does not respond to the brakes.

2) Clustering-based bad runway condition detection: The K-medoids clustering method detailed in Section III-B is performed on the Radar tracks data. Two kind of metrics are used: both SRV metrics [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF], [START_REF] Brigant | Computing distances and geodesics between manifoldvalued curves in the SRV framework[END_REF], i.e. without and with alignment, which are not supposed to take into account the 1 http://www.eurocontrol.int/articles/advanced-surface-movement-guidanceand-control-systems-smgcs deceleration of the aircraft and the slippery metric defined in Section II. Figures 9 and 10 show that using both SRV metrics, the trajectories are well clustered by the algorithm relatively to their morphologies. Indeed, cluster 2 contains mostly trajectories with a large lateral deviation, while cluster 3 corresponds to aircraft trajectories following the runway centreline (reference trajectory). Cluster 4 contains mostly standard trajectories even if there is several trajectories with slight lateral deviations. Cluster 1 corresponds to outliers trajectories. Note that, as expected, the SRV metric with alignment is more efficient for clustering trajectories morphologies than the SRV metric without alignment. We have previously seen that the latter metric may detect differences in parametrization that could influence the clustering. On the other hand, these differences in parametrization may be hidden by the shape differences and produce false alarm situations. Moreover, as explained before, if lateral deviations in landing trajectories correspond to abnormal behavior, it does not necessarily represent skidding. In practice, the detection of bad runway conditions has to take into account the deceleration profile (see Section II). It will be interesting to compare these results with the one obtained by using the slippery metric that is built for taking into account the second-order information. When adding the information of the deceleration in the metric, the results are slightly different. Indeed, from a geometric point of view, results are quite similar to the one described with the SRVF metric. However, if cluster 2 contains the most deviant trajectories, we observed one deviant trajectory in cluster 4 (Figure 11).

In-depth examination shows that it cannot be considered as a slipped trajectory. Indeed, Figure 12 compares this trajectory and its deceleration profile with the medoid of cluster 2 (deviant trajectories). Again the blue trajectory deviates from the red one but the deceleration profiles are significantly different: the deceleration level of the blue profile is higher than the red medoid deceleration profile. This means that the runway is not slippery, since there is a sufficient grip to obtain a good braking action. From an operational point of view, it seems the aircraft speed was high when arriving on exiting the runway which led on high braking. Briefly speaking, unlike the SRV metric which is only based on the geometric shape, the benefits of the metric with deceleration is twofold: firstly, it detects lateral deviations, and secondly, it removes false positives by taking into account the deceleration part. Finally, an expert study shows that there was 3.7% of slippery alerts in the dataset, among which 92.3% are in cluster 2 (red trajectories). This results confirms the good performance of the proposed metric.

V. CONCLUSION AND FUTURE WORK

The ability to detect bad runway conditions without resorting to intrusive on-site measurements will be of major interest both for economic considerations and for passenger experience. In the future, aircraft will be able to downlink their flight and taxi parameters, but until then, only radar tracks may be used. Clustering landing trajectories into skating and normal ones falls within the frame of functional data. However, some extra information must be added to really distinguish between trajectories experiencing slipping and those resulting from late breaking action of the pilot. In this work, a new distance between curves especially tailored for bad runway adherence conditions has been introduced. Using the slip angle as a measure of skating, a deformation energy is built, which associates the square derivative of the slip angle to an admissible variation of a trajectory. Using the fact that all curves are planar, its expression can be reduced to a simple semi-metric. Finally, its integration along a smooth homotopy between two curves gives an energy. Minimizing it among all possible such homotopies allows the computation of a geodesic path, with associated minimal distance. Using this new metric within a clustering algorithm yields a new methods for clustering landing trajectories that outperforms usual shape-based metrics on both simulated and real data. The cost of distance evaluation is still compatible with offline use of the procedure, in the order of ten minutes for 15,000 pairs of trajectories. It a future implementation, parallel evaluation of the distance matrix will be used, with an expected speed-up close to the number of available compute nodes.

The metric introduced is of Finsler type. The investigation of some properties of it, like its spray, was beyond the scope of the present study which focused on practical applications. It will be for further research.

  S. Puechmorel and F. Nicol are with Université de Toulouse and Ecole Nationale de l'Aviation Civile (ENAC), 31055 Toulouse, France. B. Gregorutti and C. Andrieu are with Safety Line, 75015 Paris, France.

Fig. 2 :

 2 Fig. 2: Simulated trajectories: 284 normal trajectories (blue) and 16 skating trajectories (red).

Fig. 4 :

 4 Fig. 4: Simulated circle trajectories: 15 constant speed trajectories (blue) and 15 braking trajectories (red).

Fig. 5 :

 5 Fig. 5: Simulated deviated trajectories: 15 constant speed trajectories (blue) and 15 braking trajectories (red). Low deviations (left) and high deviations (right).

Fig. 6 :

 6 Fig.6: HCA with (a) the slippery metric, (b) the SRV metric without and with alignment. Two predicted clusters: cluster 1 (blue) and cluster 2 (red).

Fig. 7 :

 7 Fig. 7: Simulated slightly deformed deviated trajectories: 15 constant speed trajectories (blue) and 15 braking trajectories (red). Low deviations (left) and high deviations (right).

Fig. 8 :

 8 Fig.8: HCA with (a) the slippery metric, (b) the SRV metric without and with alignment. Two predicted clusters: cluster 1 (blue) and cluster 2 (red).

Fig. 9 :

 9 Fig. 9: K-Medoids clustering (SRV metric without alignment).

Fig. 10 :

 10 Fig. 10: K-Medoids clustering (SRV metric with alignment).

Fig. 11 :

 11 Fig. 11: K-Medoids clustering with the slippery metric.

Fig. 12 :

 12 Fig.12:The most deviant trajectory of cluster 4 (blue), the medoid of cluster 2 (red) and their related deceleration profiles.

TABLE I :

 I number of skating and normal trajectories.

	Metric	Cluster 1	Cluster 2
		skating normal	skating normal
	Slippery	0	284	16	0
	SRV without alignment	0	284	16	0
	SRV with alignment	0	284	16	0
	ERP	0	284	16	0
	DTW	0	100	16	184
	EDR	10	284	6	0
	LCSS	10	284	6	0
	Discrete Fréchet	0	284	16	0
	Hausdorff	0	284	16	0
	SSPD	5	188	11	96

TABLE II :

 II clustering evaluation.

	Metrics	FNR	FPR Matthews
	DTW	0	0.92	0.168
	EDR/LCSS	0.034	0	0.602
	SSPD	0.259	0.90	0.164

  and Table IV . TABLE III: number of skating and normal trajectories.

	Metric		Cluster 1	Cluster 2
		skating normal	skating normal
	Slippery	0	284	16	0
	SRV without alignment	0	284	16	0
	SRV with alignment	0	284	16	0
	ERP	1	145	15	139
	Discrete Fréchet	1	145	15	139
	Hausdorff	1	145	15	139

TABLE IV :

 IV clustering evaluation for the ERP, Hausdorff and Discret Fréchet distances.

	FNR	FPR Matthews
	0.007	0.96	0.21

TABLE V :

 V number of skating and normal trajectories.

	Metric		Cluster 1	Cluster 2
		skating normal	skating normal
	Slippery	0	284	16	0
	SRV without alignment	2	186	14	198
	SRV with alignment	3	213	13	71

TABLE VI :

 VI clustering evaluation.

			Metrics				FNR	FPR Matthews
			SRV without alignment	0.011	0.93	0.141
			SRV with alignment		0.014	0.84	0.315
		0.30						0.30
		0.25						0.25
		0.20						0.20
	Y	0.15					Y	0.15
		0.10						0.10
		0.05						0.05
		0.00						0.00
		0.97	0.98	0.99	1.00	1.01		0.97	0.98	0.99	1.00	1.01
				X						X
							(a)
		0.30						0.30
		0.25						0.25
		0.20						0.20
	Y	0.15					Y	0.15
		0.10						0.10
		0.05						0.05
		0.00						0.00
		0.97	0.98	0.99	1.00	1.01		0.97	0.98	0.99	1.00	1.01
				X						X
							(b)
		0.30						0.30
		0.25						0.25
		0.20						0.20
	Y	0.15					Y	0.15
		0.10						0.10
		0.05						0.05
		0.00						0.00
		0.97	0.98	0.99	1.00	1.01		0.97	0.98	0.99	1.00	1.01
				X						X

TABLE VII :

 VII number of braking and skating curves for lowmedium deviations.

	Metric	Cluster 1	Cluster 2
		skating braking skating braking
	Slippery	0	15	15	0
	SRV without alignment	0	15	15	0
	SRV with alignment	10	10	5	5

TABLE VIII :

 VIII number of braking and skating curves for medium-high deviations.

	Metric		Cluster 1	Cluster 2
		skating braking skating braking
	Slippery	0	15	15	0
	SRV without alignment

TABLE IX :

 IX Low-medium deviations: number of braking and skating curves for slightly deformed curves.

	Metric		Cluster 1	Cluster 2
		skating braking skating braking
	Slippery	0	15	15	0
	SRV without alignment	0	11	15	4
	SRV with alignment	7	8	8	7

TABLE X :

 X clustering evaluation.

	Metrics	FNR FPR Matthews
	SRV without alignment	0	0.2	0.76
	SRV with alignment	0.47	0.47	0.07

TABLE XI :

 XI Medium-high deviations: number of braking and skating curves for slightly deformed curves.

	Metric		Cluster 1	Cluster 2
		skating braking skating braking
	Slippery	0	15	15	0
	SRV without alignment


Assembly-based Metagenomics

FM database

New species discovery

Assembly-based Metagenomics The Tchepalo strains display particular features:

-Smallest genome (<1.9 Mb) -Lack of cellobiose genes (and still uncharaterized genes) -No new functions (genes) except prophages and surface polysaccharides

Metagenomic bining at the nucleotide level

Average coverage 1/2 coverage 0 coverage

~50% variable nucleotides

This analysis shows that -at least 2 strains of L. delbrueckii are present in this sample -Their genome ~0.8% nt variability / reference genome

Assembly-based analysis

Characterization of large contigs encoding -glycosyl transferases, -surface polysaccharides, -phage genes, -mobile elements -Unknown functions

MetaGeno-Species

 Characterization of specific genes in the strains present in Tchepalo

Lactobacillus delbrueckii subsp. Jakobsenii in Tchepalo

Is Lactobacillus delbrueckii subsp. jakobsenii evolving toward Tchepalo specialisation and symbiosis?

Genomic features: has lost (only?) cellobiose sugar assimilation gene Has lost genes of unknown functions (!!!) differs by mobile elements and hypothetical proteins; No obvious gene acquisition, e.g. to allow it to adapt to sorghum fermentation (as far as our analysis shows).

Fastidious growth Always associated to Lactobacillus fermentum

Conclusion

Traditional fermented food offers new ecosystems to explore

Metagenomics is a straightforward and cheap approach to provide a real insight in their study New species and potentially new and original models of interaction / symbiosis may be characterized

The study of number of these products remains to be done in third countries and may offer them opportunities for development

Assembly Marker profiling MetaGeno- Species Lactobacillus delbrueckii subsp. ja

		~10-15 million reads
		~ 10-200 000 contigs
		55.86%
	Lactobacillus fermentum	40.22%
	Lactobacillus pontis	1.15%
	Lactobacillus agilis	0.82%
	Lactobacillus reuteri	0.74%
	Pediococcus acidilactici	0.27%
	Weissella confusa	0.22%
	Weissella cibaria	0.21%
	Saccharomyces cerevisiae	0.09%
	Arcobacter butzleri RM4018 RM40 0.03%

contigs Assembly Marker profiling MetaGeno- Species Lactobacillus delbrueckii subsp. ja

		55.86%
	Lactobacillus fermentum	40.22%
	Lactobacillus pontis	1.15%
	Lactobacillus agilis	0.82%
	Lactobacillus reuteri	0.74%
	Pediococcus acidilactici	0.27%
	Weissella confusa	0.22%
	Weissella cibaria	0.21%
	Saccharomyces cerevisiae	0.09%

~10-15 million reads ~ 10-200 000

Arcobacter butzleri RM4018 RM40 0.03% FM database New species discovery Driven culturomics Development of a Sorghum based beverage derived from Tchepalo (African beer) processing LAB fermentation Yeast fermentation No fermentation Cultural versus Metagenomic view of Tchepalo lactic step

		91%
			Lactobacillus
			fermentum
			Lactobacillus
	39%	56%	delbrueckii

2%

Cultural Metagenomics

Metagenomics indicates that beside L. fermentum, which was believed to be the major actor of the LAB fermentation, an L. delbruekii relatives may also play a major role Cultural

versus Metagenomic view of Tchepalo lactic step

			Food-Microb	16S
	Cultural	L. fermentum	56.12%	15%
	2% L. delbruekii	38.82%	82%
		L. pontis	1.43%	0%
		91% B coagulans	0.72%	0.03%
		L. agilis L. reuteri W. confusa	0.93% Lactobacillus 0.32% 1.03% 0.09% fermentum 0.24% 0.02%
	Metagenomics	S. cerevisiae	0.14%	nd
	39%	P. acidilactici W. cibaria Cyberlindnera fabianii 56% Sphingomonas melonis	0.29% Lactobacillus -0.21% -delbrueckii 0.00% nd 0.00% 0.57%

-16S and shotgun metagenomics give slightly different results but:  L. delbruekii confirmed as a major actor in Tchepalo processing although underscored by cultural approach Metagenomic

profiling of Tchepalo L. delbrueckii strains

	Lactobacillus delbrueckii	Lactobacillus delbrueckii
	subsp. lactis KCCM 34717	subsp. jakobsenii DSM 26046

-Metagenomic profiling allowed to point out the L. delbruekii subspecies genome, here a strain isolated from a comparable product in a neighbor country -Modification of growth condition allowed isolating definitivly strains from Tchepalo Comparative genomics of L. delbrueckii strains subsp. lactis subsp. delbruekii subsp. bulgaricus subsp. indicus subsp. jakosenii

Lactobacillus delbrueckii subsp. lactis KCCM 34717Lactobacillus delbrueckii subsp. jakobsenii DSM 26046


Imad Ammoun

Caroline I Kothe

Catherine Béal

Rita Yaacoub

C Beal

P Renault

Metagenomics reveal the specificities of the microbiota of three traditional Lebanese goat milk fermented products

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Conclusions

Product selection, production and producer's survey

Results and discussion

Materials and Methods

In Lebanon, the production of traditional fermented goat dairy products is carried out at small scale, and only in a limited number of production units. As a result, there is a real risk of extinction of these products and loss of corresponding microbial resources (Dimassi et al 2020). Therefore, the investigation of their microbial system is becoming urgent matter. As a response, three products were selected for this study, namely Ambriss, Serdaleh and Labneh El Darff. This study aims to reveal their characteristics and microbial diversity in order to preserve their biological resources for further evaluation for their technological properties. According to the survey, the production is very artisanal, mainly depending on women housekeepers (60%) aged between 40 and 78 years old and displays limited access to local market place (38%).

The examination of culture results showed an important variability in the different enumerations. A dominance of mesophilic lactobacilli (5.4 -6.8 Log CFU/g) over lactococci (0.5 -5.2 Log CFU/g) was observed in the three products. Thermophilic counts (MRS at 42 °C) were highly variable in Labneh El Darff (1 -6.2 Log CFU/g) compared to its other counterparts. Plate counts of lactobacilli at 30 °C were more pronounced (P < 0.05) in Labneh El Darff (6.8 Log CFU/g) than those recorded for Ambriss and Serdaleh (5.4 and 6.0 Log CFU/g respectively).

Species determination by 16S sequencing of the isolates showed that L. plantarum was a dominant species in Ambriss (38%) followed by L. rhamnosus and Enterococcus durans (18% for each). Serdaleh presents 10 different species similar to Ambriss, in which L. plantarum was predominant (63%), whereas in Labneh El Darff, Lb. brevis was the most abundant followed by L. bulgaricus subsp. bulgaricus, and Lactococcus lactis (16% for each).

Introduction

Culture-independent approaches

Producer survey Amphora Goat skin Using both culture-dependent and independent approaches, we demonstrated that Serdaleh and Ambriss possessed an original living microbiota that was present to a lesser extent in Labneh El Darff. The first two products share an association of Lactococcus lactis and three Lactobacillus species, with different proportions as they originate from 2 distinct regions, even if they share similar ways of processing. The detection of contaminating bacteria in some products focuses the attention on the implementation of Good Manufacturing Practices and Good Hygienic Practices in such artisanal productions. In addition, a commercialization policy is needed to develop other social blocks that enable housekeepers to be involved.

The estimation of microbial complexity by metagenetic study, based on the 4 indices of alfa diversity (Figure 2) showed that Labneh El Darff presents a low bacterial diversity (P < 0.05) and is dominated by Lactobacillus delbrueckii subsp. bulgaricus, originating from Laban.

Culture-dependent approaches

  Serhan M. and Mattar . 2013. Characterization of four artisanal goat milk cheese: Darfieh, Ariche, Shanklish and Serdal by physico-chemical, microbiological and sensory analyses. Journal of Food Agricultureand Environment Vol 11 (3 and 4): 97 -101. Dimassi O., Iskandarani Y., Afram M., Akiki R., Rached M. 2020. Production and physicochemical properties of Labneh anbaris, a traditional fermented cheese like product in Lebanon. International Journal of Environment, Agriculture and Biotechnology Vol 5 (3): 509 -516. Traditional Lebanese cheese in danger of disappearing https://www.beirut.com/l/56124, 17/08/2018 SRVA, 2005. Proposal for the protection of geographical indication in Lebanon. Report on the intention of IDEAS centre Arabesque10, 1204 Geneva. 33p.

  on different media while their DNA was extracted for the 16S metagenetic and Shotgun metagenomic analyses. From 380 isolates, 279 strains were genetically identified by 16S rDNA. Metagenomic microbial composition was estimated (i) by mapping the reads against the representative clade-specific marker catalog using the MetaPhlAn tool and (ii) by using YchF gene marker for binning of the assembled metagenomic contigs.. A questionnaire regarding farming and production practices was conducted on 50 producers, from which 15 samples of Serdaleh, 11 of Ambriss and 4 of Labneh El Darff were collected across three different regions.

Figure 1 .

 1 Figure 1. Schematic flow chart for the three products related to their geographical origin on the Lebanese map.

Figure 3 :

 3 Figure 3 : Taxonomic composition plot of the relative abundance of the 20 most abundant bacterial OUTs assigned at the species level in Ambriss, Serdaleh and Labneh El Darff products.

Table 1 :

 1 Surveys

	Ambriss Serdaleh	Labneh El Darff

' results related to production and farming practices of the three products Exploring beta diversity, Serdaleh and Ambriss presented 2 major clusters sharing a similar microbiota mainly composed of Lactobacilli. . Combined analysis with metagenomics showed that Serdaleh has a distinct homogenous clade dominated by Lactococcus lactis subsp lactis (> 65% of reads) followed by the species L. kefiranofaciens, L. kefiri and L. diolivorans. Sample A-4 was dominated by Streptococcus parasuis group while sample A-6 showed high amount of L. gasseri (~36%), L. rhamnosus (~23%) and Enterobacteriaceae (Figures

2 and 3

). All Serdaleh samples were contaminated by Mycoplasmopsis agalactiae and a low level of Listeria monocytogenes was detected in 4 Ambriss and 2 Serdaleh. Capra hircus reads, corresponding to goat DNA were shown in abundance in all Labneh samples, probably originating from goat's skin during processing.


I Ammoun

C I Kothe

C Beal

R Yaacoub

P Renault

Revealing the microbial heritage of traditional Lebanese goat milk fermented products

In Lebanon, the production of traditional fermented goat dairy products is carried out at small scale, and only in a limited number of production units. As a result, there is a real risk of extinction of these products and loss of corresponding microbial resources. The understanding of their microbial ecosystems is becoming necessary and urgent. The aim of this work was to provide an overview of their microbial biodiversity, to isolate and characterize representative strains, and to point out specificities of these microbiota by combining approaches, including assembly-and deep profilingmetagenomics.

We will highlight here the results of "Ambriss" and "Serdaleh" produced respectively in the Bekaa and Chouf regions, by successively filling an amphora with raw goat's milk, draining the whey and salting the surface until it is filled. The production is found to be typically artisanal, mainly depending on women housekeepers (60%) aged between 40 and 78 years old and displays limited access to local market place (38%).

Producer's survey

Microbial composition

Figure 4 : Composition based on OTUs Culture-dependent and -independent analyses gave drastically different results. Metagenomics showed that Serdaleh and Ambriss share a similar microbiota but with different proportions of L. lactis and L. kefiranofaciens, L. kefiri and L. diolivorans. Differences in local conditions (temperature, salt…) may explain the shift in population in these products. This composition is similar to that of kefir, a probiotic product, whose microbiota is driven by a particular structure in polysaccharide aggregates, the kefir grain.

In addition, all Serdaleh samples contained Mycoplasmopsis agalactiae, the causing agent of contagious agalactia in sheep and goats.

CONCLUSION & PERSPECTIVES

Ambris and Serdaleh, emblematic products of several millennia and considered part of the Lebanese heritage, share the same microbial composition as kefir, a popular probiotic product worldwide. This fact adds to the interest of preserving their production. However, there are several issues that need to be resolved before further development of these products can take place. ABR genes were detected in most samples, particularly in Serdaleh, whose dominant L. lactis strains contain a plasmid with a multiresistance island.

-Additional analyses are required to understand the resilience of the ecosystem established in the amphora, which are used only during goat milk production season.

-Sanitary practices have to be improved for milk production -Finally, appropriate social and marketing policies should be developed and implemented to support and maintain the main actors of these productions RESULTS

Fig. 2 .Fig. 3 .

 23 Fig. 2. Overview of methodological approaches used in this study.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Example of phylogenetic analysis of product MAGs and reference strain genomes.

Exudation Laban Salt Milk Salt Drainage Labneh El Darf Ambriss Serdaleh Periodic filling Amphora Goat skin Fig 1.

  Schematic flow chart for three emblematic products related to their geographical origin

					Ambriss	Serdaleh	Labneh El Darff
				Producers dedicated to only to 1 dairy product (%)	80	51.61	100
				Extensive goat Farming system (%)	50	66.70	66
				Milk processed by farmers-producers (L/year)	160-3200	220-7200	135-1350
				Milk processed by producers only (L/year)	140-1890	110-1250	135-850
				Average yearly cheese production (Kg)	108.8	411.9	308.3
				Production yield ( Kg/100L milk)	19.44	18.89	16.35
				Producers (%) implementing a processing period of 1-3 months	0	15	55.6
				Producers (%) implementing a processing period of 4-6 months	100	75	44.4
			WORKPLAN
	A survey regarding farming and production practices
	was conducted on 50 producers, from which 15
	samples of Serdaleh, 11 of Ambriss and 4 of Labneh
	El Darff were collected across the three different
	regions.
	Samples were cultured on different media and 279
	strains were genetically identified by 16S rDNA.
	Cheese DNA was extracted for the 16S metagenetic
	and	Shotgun	metagenomic	sequencing.
	Bioinformatics was performed as described below.

Table 1 :

 1 Survey's results related to production and farming practices

Conseils de mise en page