L Roques

M Cristofol

On the determination of the nonlinearity from localized measurements in a reaction-diffusion equation

Keywords: reaction-diffusion, heterogeneous media, uniqueness, inverse problem

This paper is devoted to the analysis of some uniqueness properties of a classical reaction-diffusion equation of Fisher-KPP type, coming from population dynamics in heterogeneous environments. We work in a one-dimensional interval (a, b) and we assume a nonlinear term of the form u (µ(x)γu) where µ belongs to a fixed subset of C 0 ([a, b]). We prove that the knowledge of u at t = 0 and of u, u x at a single point x 0 and for small times t ∈ (0, ε) is sufficient to completely determine the couple (u(t, x), µ(x)) provided γ is known. Additionally, if u xx (t, x 0 ) is also measured for t ∈ (0, ε), the triplet (u(t, x), µ(x), γ) is also completely determined. Those analytical results are completed with numerical simulations which show that, in practice, measurements of u and u x at a single point x 0 (and for t ∈ (0, ε)) are sufficient to obtain a good approximation of the coefficient µ(x). These numerical simulations also show that the measurement of the derivative u x is essential in order to accurately determine µ(x).

Introduction and ecological background

Reaction-diffusion models (hereafter RD models), although they sometimes bear on simplistic assumptions such as infinite velocity assumption and completely random motion of animals [START_REF] Holmes | Are diffusion models too simple? A comparison with telegraph models of invasion[END_REF], are not in disagreement with certain dispersal properties of populations observed in natural as well as experimental ecological systems, at least qualitatively [START_REF] Shigesada | Biological invasions: theory and practice[END_REF][START_REF] Turchin | Quantitative analysis of movement: measuring and modeling population redistribution in animals and plants[END_REF][START_REF] Murray | Mathematical Biology[END_REF][START_REF] Okubo | Diffusion and ecological problems -modern perspectives[END_REF]. In fact, since the work of Skellam [START_REF] Skellam | Random dispersal in theoretical populations[END_REF], RD theory has been the main analytical framework to study spatial spread of biological organisms, partly because it benefits from a well-developed mathematical theory.

The idea of modeling population dynamics with such models has emerged at the beginning of the 20 th century, with random walk theories of organisms, introduced by Pearson and Blakeman [START_REF] Pearson | Mathematical contributions of the theory of evolution. A mathematical theory of random migration[END_REF]. Then, Fisher [START_REF] Fisher | The wave of advance of advantageous genes[END_REF] and Kolmogorov, Petrovsky, Piskunov [START_REF] Kolmogorov | Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF] independently used a reaction-diffusion equation as a model for population genetics. The corresponding equation is

∂u ∂t -D ∂ 2 u ∂x 2 = u (µ -γu), t > 0, x ∈ (a, b) ⊂ R, (1.1)
where u = u(t, x) is the population density at time t and space position x, D is the diffusion coefficient, and µ and γ respectively correspond to the constant intrinsic growth rate and intraspecific competition coefficients. In the 80's, this model has been extended to heterogeneous environments by Shigesada et al. [START_REF] Shigesada | Traveling periodicwaves in heterogeneous environments[END_REF]. The corresponding model is of the type:

∂u ∂t - ∂ ∂x D(x) ∂u ∂x = u (µ(x) -γ(x)u), t > 0, x ∈ (a, b). (1.2)
The coefficients µ(x) and γ(x) now depend on the space variable x and can therefore include some effects of environmental heterogeneity. More recently, this model revealed that the heterogeneous character of the environment played an essential role on species persistence and spreading, in the sense that for different spatial configurations of the environment, a population can survive or become extinct and spread at different speeds, depending on the habitat spatial structure ( [START_REF] Shigesada | Biological invasions: theory and practice[END_REF], [START_REF] Berestycki | Analysis of the periodically fragmented environment model: I -Species persistence[END_REF], [START_REF] Cantrell | Spatial ecology via reaction-diffusion equations[END_REF], [START_REF] El Smaily | Homogenization and influence of fragmentation in a biological invasion model[END_REF], [START_REF] Hamel | Spreading speeds in slowly oscillating environments[END_REF] , [START_REF] Roques | Mathematical analysis of the optimal habitat configurations for species persistence[END_REF], [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF]). Thus, determining the coefficients in model (1.2) is an important question, even for areas other than ecology (see [START_REF] Xin | Front propagation in heterogeneous media[END_REF] and references therein).

In this paper, we focus on the case of constant coefficients D and γ:

∂u ∂t -D ∂ 2 u ∂x 2 = u (µ(x) -γu), t > 0, x ∈ (a, b), (1.3)
and we address the question of the uniqueness of couples (u, µ(x)) and triples (u, µ(x), γ) satisfying (1.2), given a localized measurement of u. Uniqueness results of this type have been obtained for reaction-diffusion models, through the Lipschtiz stability of the coefficient with respect to the solution u. Lipschtiz stability is generally obtained by using the method of Carleman estimates [START_REF] Bukhgeim | Global uniqueness of a class of multidimensional inverse problems[END_REF]. Several publications starting from the paper by Isakov [START_REF] Isakov | An uniqueness in inverse problems for semilinear parabolic equations[END_REF] and including the recent overview of the method of Carleman estimates applied to inverse coefficients problems [START_REF] Klibanov | Carleman estimates for coefficient inverse problems and numerical applications[END_REF] provide results for the case of multiple measurements. The particular problem of the uniqueness of the couple (u, µ(x)) satisfying (1.3) given such multiple measurements has been investigated, together with Lipschtiz stability, in a previous work [START_REF] Cristofol | Biological invasions: Deriving the regions at risk from partial measurements[END_REF]. Placing ourselves in a bounded domain Ω of R N with Dirichlet boundary conditions, we had to use the following measurements: (i) the density u(0, x) in Ω at t = 0; (ii) the density u(t, x) for (t, x) ∈ (t 0 , t 1 ) × ω, for some times 0 < t 0 < t 1 and a subset ω ⊂⊂ Ω; (iii) the density u(θ, x) for all x ∈ Ω, at some time θ ∈ (t 0 , t 1 ).

Although the result of [START_REF] Cristofol | Biological invasions: Deriving the regions at risk from partial measurements[END_REF] allows to determine µ(x) using partial measurements of u(t, x), assumption (iii) implies that u has to be known in the whole set Ω. This last measurement (iii) is a key assumption in several other papers on uniqueness and stability of solutions to parabolic equations with respect to parameters (see Imanuvilov and Yamamoto [START_REF] O Y Immanuvilov | Lipschitz stability in inverse parabolic problems by the Carleman estimate[END_REF], Yamamoto and Zou [START_REF] Yamamoto | Simultaneous reconstruction of the initial temperature and heat radiative coefficient[END_REF], Belassoued and Yamamoto [START_REF] Belassoued | Inverse source problem for a transmission problem for a parabolic equation[END_REF] for scalar equations and Cristofol, Gaitan and Ramoul [START_REF] Cristofol | Inverse problems for a two by two reaction-diffusion system using a carleman estimate with one observation[END_REF] or Benabdallah, Cristofol, Gaitan and Yamamoto [START_REF] Benabdallah | Inverse problem for a parabolic system with two components by measurements of one component[END_REF] for systems).

Here, contrarily to previous results obtained for this type of reactiondiffusion models, there are some regions in (a, b) where u is never measured: we only require to know (i') the density u(0, x) in (a, b) at t = 0 and (ii') the density u(t, x 0 ) and its spatial derivative ∂u ∂x (t, x 0 ) for t ∈ (0, ε) and some point x 0 in (a, b) (see Remark 2.5 for a particular example of hypothesis (ii')). Thus a measurement of type (iii) is no more necessary. Furthermore, we show simultaneous uniqueness of two coefficients µ(x) and γ provided that measurements of the second derivative ∂ 2 u ∂x 2 (t, x 0 ) are available. Our paper is organized as follows: in the next section, we give precise statements of our hypotheses and results; Section 3 is then dedicated to the proof of the results. Section 4 is devoted to the description of numerical examples illustrating how the coefficient µ(x) can be approached using measures of the type (i') and (ii'). Those results are further discussed in Section 5.

Hypotheses and main results

Let (a, b) be an interval in R. We consider the problem:

       ∂u ∂t -D ∂ 2 u ∂x 2 = u (µ(x) -γu), t ≥ 0, x ∈ (a, b), α 1 u(t, a) -β 1 ∂u ∂x (t, a) = 0, t > 0, α 2 u(t, b) + β 2 ∂u ∂x (t, b) = 0, t > 0, u(0, x) = u i (x), x ∈ (a, b). (P µ,γ )
Our hypotheses on the coefficients are the following. Firstly, we assume that:

µ ∈ M := {ψ ∈ C 0,η ([a, b]) such that ψ is piecewise analytic on (a, b)}, (2.4
) for some η ∈ (0, 1]. The space C 0,η corresponds to Hölder continuous functions with exponent η (see e.g. [START_REF] Friedman | Partial differential equations of parabolic type[END_REF]). A function ψ ∈ C 0,η ([a, b]) is called piecewise analytic if it exists n > 0 and an increasing sequence (i k ) 1≤k≤n such that i 1 = a, i n = b, and

for all x ∈ (a, b), ψ(x) = n-1 j=1 χ [i j ,i j+1 ) (x)ϕ j (x),
for some analytic functions ϕ j , defined on the intervals [i j , i j+1 ], and where χ [i j ,i j+1 ) are the characteristic functions of the intervals [i j , i j+1 ) for j = 1, . . . , n -1.

We also assume that γ is a positive constant and that the boundary coefficients satisfy:

α 1 , α 2 , β 1 , β 2 ≥ 0 with α 1 + β 1 > 0 and α 2 + β 2 > 0.
(2.5)

We furthermore make the following hypotheses on the initial condition:

u i ≥ 0, u i ≡ 0 and u i ∈ C 2,η ([a, b]), (2.6)
for some η in (0, 1), that is u i is a C 2 function such that u ′′ i is Hölder continuous. In addition to that, we assume the following compatibility conditions:

α 1 u i (a) -β 1 u ′ i (a) = 0, α 2 u i (b) + β 2 u ′ i (b) = 0, δ β 1 u ′′ i (a) = 0, δ β 2 u ′′ i (b) = 0, (2.7
) where δ y is verifies: δ 0 = 1 and δ y = 0 if y = 0. We also need to assume that:

measure({x ∈ (a, b), u i (x) = 0}) = 0. (2.8)
Under the assumptions (2.4)-(2.7), for each µ ∈ M and γ > 0, the problem (P µ,γ ) has a unique solution u ∈ C 2,η 1,η/2 ([0, +∞) × [a, b]) (i.e. the derivatives up to order two in x and order one in t are Hölder continuous, see [START_REF] Friedman | Partial differential equations of parabolic type[END_REF][START_REF] Pao | Nonlinear Parabolic and Elliptic Equations[END_REF] for a definition of Hölder continuity). Existence, uniqueness and regularity of the solution u are classical. See e.g. [START_REF] Pao | Nonlinear Parabolic and Elliptic Equations[END_REF]Ch. 1].

Let us state our main results: Theorem 2.1. Let µ, μ ∈ M, and γ > 0, and assume that the solutions u and ũ to (P µ,γ ) and (P μ,γ ) satisfy, at some x 0 ∈ (a, b), and for some ε > 0 and all t in (0, ε):

u(t, x 0 ) = ũ(t, x 0 ), (2.9) ∂u ∂x (t, x 0 ) = ∂ ũ ∂x (t, x 0 ). ( 2

.10)

Assume furthermore that

u i (x 0 ) = 0 or ∂ 2 u ∂x 2 (t, x 0 ) = ∂ 2 ũ ∂x 2 (t, x 0 ) for t ∈ (0, ε). (2.11) Then, we have µ ≡ μ on [a, b] and consequently u ≡ ũ in [0, +∞) × [a, b]. If β 1 > 0 (resp. β 2 > 0), this statement remains true when x 0 = a (resp. x 0 = b). Remark 2.2. This result remains valid if γ = γ(x) is a given, positive function in C 0,η ([a, b]).
However, the conclusion of Theorem 2.1 is not true in general without the assumption (2.10): 2 ) = ũ(t, a+b 2 ) for all t ≥ 0. Under an additional assumption on the initial condition u i , we are able to obtain a uniqueness result for triples (u, µ, γ): Theorem 2.4. Let µ, μ ∈ M, and γ, γ > 0. Assume that, at some x 0 ∈ (a, b), u i (x 0 ) = 0. Assume furthermore that the solutions u and ũ to (P µ,γ ) and (P μ,γ ) satisfy, for some ε > 0 and for all t in (0, ε): Remarks 2.5.

u(t, x 0 ) = ũ(t, x 0 ), (2.12) ∂u ∂x (t, x 0 ) = ∂ ũ ∂x (t, x 0 ), (2.13) ∂ 2 u ∂x 2 (t, x 0 ) = ∂ 2 ũ ∂x 2 (t, x 0 ). ( 2
• A particular example where hypotheses (2.9-2.11) of Theorem 2.1 (resp. hypotheses (2.12-2.14) of Theorem 2.4) are fulfilled is whenever, for some subset ω of (a, b), u(t, x) = ũ(t, x) for t ∈ (0, ε) and all x ∈ ω (resp. x 0 ∈ ω and u(t, x) = ũ(t, x) in (0, ε) × ω). Note that, under this hypothesis, the previous results [START_REF] Cristofol | Biological invasions: Deriving the regions at risk from partial measurements[END_REF] did not imply uniqueness; indeed, an additional assumption of type (iii) was required (cf. the introduction section).

• The uniqueness result of Theorem 2.4 cannot be adapted to the stationary equation associated to (P µ,γ ): -p ′′ = p (µ(x)γp) (see e.g. [START_REF] Berestycki | Analysis of the periodically fragmented environment model: I -Species persistence[END_REF] for the existence and uniqueness of the stationary state p > 0). Indeed, for any τ ∈ (0, 1), setting μ = µτ γp and γ = (1τ )γ, we obtain -p ′′ = p (μ(x)γp), whereas μ ≡ µ and γ ≡ γ. Thus, a measurement of p, even on the whole interval [a, b], does not provide a unique couple (µ, γ).

•

Proofs

Let µ, μ ∈ M, and γ, γ > 0. Let u be the solution to (P µ,γ ) and ũ the solution to (P μ,γ ). We set

U := u -ũ and m := µ -μ.
The function U satisfies:

∂U ∂t -D ∂ 2 U ∂x 2 = μU -γU (u + ũ) + u(m -u(γ -γ)), (3.15)
for t ≥ 0 and x ∈ (a, b), and (3.17)

α 1 U (t, a) -β 1 ∂U ∂x (t, a) = 0, α 2 U (t, b) + β 2 ∂U ∂x (t, b) = 0, t > 0, U (0, x) = 0, x ∈ (a, b). ( 3
Step 1: We prove that m(x 0 ) = 0.

It follows from hypothesis (2.9) that, for all t ∈ [0, ε), U (t, x 0 ) = 0 and thereby, ∂U ∂t (t, x 0 ) = 0 for all t ∈ [0, ε).

If u i (x 0 ) = 0, then, since U (0, •) ≡ 0 we deduce from (3.17) applied at t = 0 and x = x 0 that u i (x 0 )m(x 0 ) = 0, and therefore m(x 0 ) = 0. If u i (x 0 ) = 0, from (2.11), we have ∂ 2 U ∂x 2 (t, x 0 ) = 0 for all t ∈ [0, ε). Applying equation (3.17) at t = ε/2 and x = x 0 , we get

u ε 2 , x 0 m(x 0 ) = 0. If x 0 ∈ (a, b
), the strong parabolic maximum principle (Corollary 5.2) applied to u implies that u(ε/2, x 0 ) > 0. As a consequence we again get m(x 0 ) = 0. Lastly, if x 0 = a and β 1 > 0 the Hopf's Lemma applied to u again implies that u(ε/2, x 0 ) > 0. Indeed, assume on the contrary that u(ε/2, x 0 ) = u(ε/2, a) = 0. The boundary condition α 1 u(t, a)β 1 ∂u ∂x (t, a) = 0 implies:

β 1 ∂u ∂x ε 2 , a = 0,
which is impossible from Hopf's Lemma (Corollary 5.2 and Theorem 5.1 (b) and (c)). Thus u(ε/2, x 0 ) > 0 and, again, m(x 0 ) = 0. A similar argument holds for x 0 = b, whenever β 2 > 0.

Under the assumptions of Theorem 2.1, we therefore always obtain m(x 0 ) = 0.

Step 2: We prove that m ≡ 0.

Let us now set

b 1 := sup{x ∈ [x 0 , b] s.t. m has a constant sign on [x 0 , x]}.
By "constant sign" we mean that either m ≥ 0 on [x 0 , x] or m ≤ 0 on [x 0 , x]. Then, four possibilities may arise: Since µ and μ belong to M , the function m also belongs to M and is therefore piecewise analytic on (a, b). Thus, the set {x ∈ (a, b) s.t. m(x) = 0} has a finite number of connected components. This contradicts (3.18) and rules out possibility (i). Now assume (ii). By continuity of m, and from hypothesis (2.8) on u i , we can assume that u i (x 1 ) > 0. Since m(x 1 ) > 0 and U (0, •) ≡ 0, it follows from (3.17) that ∂U ∂t (0, x 1 ) = u i (x 1 )m(x 1 ) > 0.

• (i) m = 0 on [x 0 , b 1 ] and b 1 < b, • (ii) m ≥ 0 on [x 0 , b 1 ] and it exists x 1 ∈ (x 0 , b 1 ) such that m(x 1 ) > 0, • (iii) m ≤ 0 on [x 0 , b 1 ] and it exists x 1 ∈ (x 0 , b 1 ) such that m(x 1 ) < 0, • (
Thus, for ε 1 > 0 small enough, U (t, x 1 ) > 0 for t ∈ (0, ε 1 ]. As a consequence, U satisfies:

   ∂U ∂t -D ∂ 2 U ∂x 2 -(μ -γu -γ ũ)U ≥ 0, t ∈ (0, ε 1 ], x ∈ (x 0 , x 1 ), U (t, x 0 ) = 0 and U (t, x 1 ) > 0, t ∈ (0, ε 1 ], U (0, x) = 0, x ∈ (x 0 , x 1 ). (3.19)
Moreover, Lemma 3.1. We have U (t, x) > 0 in (0, ε 1 ) × (x 0 , x 1 ).

Proof of Lemma 3.1: Set W = U e -λt , for some λ > 0 large enough such that c(t, x)

:= μ -γu -γ ũ -λ ≤ 0 in (x 0 , x 1 ). The function W satisfies ∂W ∂t -D ∂ 2 W ∂x 2 -c(t, x)W ≥ 0, t ∈ (0, ε 1 ], x ∈ (x 0 , x 1 ).
Assume that it exists a point (t * , x * ) in (0, ε 1 )×(x 0 , x 1 ) such that U (t * , x * ) < 0. Then, since W (t, x 0 ) = 0 and W (t, x 1 ) > 0 for t ∈ (0, ε 1 ), and since

W (0, x) = U (0, x) = 0, W admits a minimum m * < 0 in (0, ε 1 ] × (x 0 , x 1 ). Theorem 5.1 (a) applied to W implies that W ≡ m * < 0 on [0, ε 1 ] × [x 0 , x 1 ], which is impossible. Thus W ≥ 0 in [0, ε 1 ]×[x 0 , x 1 ]
. Theorem 5.1 (a) and (c) then implies that W > 0 and consequently U (t, x) > 0 in (0, ε 1 ) × (x 0 , x 1 ).

Since U (t, x 0 ) = 0, the Hopf's lemma (Theorem 5.1 (b) and (c)) also implies that ∂U ∂x (t, x 0 ) > 0 for all t ∈ (0, ε 1 ). This contradicts hypothesis (2.10). Possibility (ii) can therefore be ruled out.

Applying the same arguments to -U , possibility (iii) can also be rejected. Finally, only (iv) remains.

Setting

a 1 := inf{x ∈ [a,
u(t, x 0 ) (m(x 0 ) -u(t, x 0 )(γ -γ)) = 0 for t ∈ [0, ε).
If x 0 ∈ (a, b), the strong parabolic maximum principle (Corollary 5.2) implies that u(t, x 0 ) > 0 for all t > 0. This remains true if x 0 = a (if β 1 > 0) or x 0 = b (if β 2 > 0); cf. the proof of Theorem 2.1. We therefore get: m(x 0 ) = u(t, x 0 )(γγ) for t ∈ (0, ε).

(3.20)

From the continuity of t → u(t, x 0 ) up to t = 0, we have m(x 0 ) = u i (x 0 )(γγ). Thus, u i (x 0 ) = 0 implies that m(x 0 ) = 0 which in turns implies from (3.20), and since u(t, x 0 ) > 0 for t > 0, that γ = γ. The end of the proof is therefore similar to that of Theorem 2.1.

Remark 3.2. Extension of the arguments used in the previous proof to higher dimensions is not straightforward. Indeed, placing ourselves in a bounded domain Ω of R N , with N ≥ 2, we may consider the largest region Ω 1 in Ω, containing x 0 and such that m has a constant sign in Ω 1 . Consider in the above proof the possibility (ii) N (instead of (ii)): m ≥ 0 on Ω 1 and it exists x 1 ∈ Ω such that m(x 1 ) > 0. Then it exists a subset ω 1 of Ω 1 , such that x 0 ∈ ∂ω 1 and u(t, x) > 0 on a portion of ∂ω 1 . However, we cannot assert that U (t, x) ≥ 0 on ∂ω 1 , and (ii) N can therefore not be ruled out as we did for (ii).

Proof of Proposition 2.3: Under the assumptions of Proposition 2.3, we observe that ũ(t, b -(xa)) is a solution of (P µ,γ ). In particular, by uniqueness, we have

u(t, x) = ũ(t, b -(x -a)), for all t ≥ 0 and x ∈ [a, b].

It follows that u(t, a+b

2 ) = ũ(t, a+b 2 ) for all t ≥ 0.

Numerical computations

The purpose of this section is to verify numerically that the measurements (2.9-2.10) of Theorem 2.1 allow to obtain a good approximation of the coefficient µ(x), when γ is known.

Assuming that µ belongs to a finite-dimensional subspace E ⊂ M and measuring the distance between the measurements of the solutions of (P µ,γ ) and (P μ,γ ) through the function

G µ (μ) = u(•, x 0 ) -ũ(•, x 0 ) L 2 (0,ε) + ∂u ∂x (•, x 0 ) - ∂ ũ ∂x (•, x 0 ) L 2 (0,ε) ,
we look for the coefficient µ(x) as a minimizer of the function G µ . Indeed, G µ (µ) = 0 and, from Theorem 2.1, this is the unique global minimum of G µ in M . Solving (P µ,γ ) by a numerical method (see Appendix B) gives an approximate solution u h . In our numerical tests, we therefore replace G µ by the discretized functional

G µ (μ) := u h (•, x 0 ) -ũh (•, x 0 ) L 2 (0,ε) + ∂u h ∂x (•, x 0 ) - ∂ ũh ∂x (•, x 0 ) L 2 (0,ε) .
Remark 4.1. Since (P µ,γ ) and (P μ,γ ) are solved with the same (deterministic) numerical method, we have G µ (µ) = 0. Thus µ is a global minimizer of G µ . However, this minimizer might not be unique.

State space E

We fix (a, b) = (0, 1) and we assume that the function µ belongs to a subspace E ⊂ M defined by:

E := μ ∈ C 0,η ([0, 1]) | ∃ (h i ) 0≤i≤n ∈ R n+1 , μ(x) = n i=0 h i • j ((n -2) (x -c i )) on [0, 1] ,
with c i = i-1 n-2 and j(x) = exp 4x 2 x 2 -4 , if x ∈ (-2, 2), 0 otherwise.

Minimization of G µ in E

For the numerical computations, we fixed D = 0.1, γ = 1, α 1 = α 2 = 0 and β 1 = β 2 = 1 (Neumann boundary conditions). Besides, we assumed that u i ≡ 0.2, ε = 0.3 and x 0 = 2/3. The integer n was set to 10 in the definition of E.

Numerical computations were carried out for 100 functions µ k in E :

µ k = n i=0 h k i • j [(n -2) (x -c i )] , k = 1 . . . 100,
whose components h k i were randomly drawn from a uniform distribution in (-5, 5).

Minimizations of the functions G µ k were performed using MATLAB's fminunc solver1 . This led to 100 functions µ * k in E, each one corresponding to a computed approximation for a minimizer of the function G µ k . In our numerical tests, we obtained values of G µ k (µ * k ) in (5 • 10 -7 , 10 -5 ), with an average of 5 • 10 -6 and a standard deviation of 2 • 10 -6 .

The values

µ k -µ * k L 2 ([0,1]) / µ k L 2 ([0,1]
) , for k = 1 . . . 100, are comprised between 5 • 10 -3 and 0.16, with an average value of 0.04 and a standard deviation of 0.03.

Fig. 1 (a) depicts an example of function µ in E, together with a function µ * which was obtained by minimizing G µ .

Test of another criterion H µ

In this section, we illustrate that measurement (2.9) alone cannot be used for reconstructing µ. Replacing G µ by:

H µ (μ) = u(•, x 0 ) -ũ(•, x 0 ) L 2 (0,ε) ,
and setting H µ (μ) := u h (•, x 0 )ũh (•, x 0 ) L 2 (0,ε) , we performed the same analysis as above, with the same samples µ k ∈ E and the same parameters.

The corresponding values of H µ (µ * k ) are comparable to those obtained in Section 4.

µ -µ * L 2 ([0,1]) / µ L 2 ([0,1]) = 0.47 and H µ (µ * ) = 2 • 10 -6 .
average 2•10 -6 , and standard deviation 3•10 -6 . However, the corresponding values of the distance µ ]) are far larger than those obtained in Section 4.2: these values are comprised between 0.08 and 1.64, with an average of 0.56 and a standard deviation of 0.34.

k -µ * k L 2 ([0,1]) / µ k L 2 ([0, 1
Using the same sample µ ∈ E as in Fig. 1 (a), we present in Fig. 1 (b) the approximation µ * obtained by minimizing H µ . In this case, the distance µµ

* L 2 ([0,1]) / µ L 2 ([0,1]) is 18 times larger than µ -µ * L 2 ([0,1]) / µ L 2 ([0,1]) .

Discussion

Studying the reaction-diffusion problem (P µ,γ ) with a nonlinear term of the type u (µ(x)γu), we have proved in Section 2 that knowing u and its first spatial derivative at a single point x 0 and for small times t ∈ (0, ε) is sufficient to completely determine the couple (u(t, x), µ(x)). Additionally, if the second spatial derivative is also measured at x 0 for t ∈ (0, ε), the triplet (u(t, x), µ(x), γ) is also completely determined. These uniqueness results are mainly the consequences of Hopf's Lemma and of an hypothesis on the set M of coefficients which µ(x) belongs to. This hypothesis implies that two coefficients in M can be equal only over a set having a finite number of connected components.

The theoretical results of Section 2 suggest that the coefficients µ(x) and γ can be numerically determined using only measurements of the solution u of (P µ,γ ) and of its spatial derivatives at one point x 0 , and for t ∈ (0, ε). Indeed, the numerical computations of Section 4 show that, when γ is known, the coefficient µ(x) can be estimated by minimizing a function G µ . The function ũ being the solution of (P μ,γ ), we defined G µ (μ) as the distance between (u, ∂u/∂x)(•, x 0 ) and (ũ, ∂ ũ/∂x)(•, x 0 ), in the L 2 (0, ε) sense.

The numerical computations presented in Section 4.2 were carried out on 100 samples of functions µ k chosen in a finite-dimensional subspace of M . In each case, a good approximation µ * k of µ k was obtained. The average relative L 2 -error between µ k and µ * k is 30 times smaller than the average relative L 2 -error between µ k and the constant function µ k (x 0 ). Thus, a measurement of u and of its first spatial derivative at a point x 0 (and for t ∈ (0, ε)) indirectly gives more information on the global shape of µ than a direct measure of µ at x 0 . These good results, in spite of the computational error, indicate L 2 -stability of the coefficient µ with respect to single-point measurements of the solution u of (P µ,γ ) and of its spatial derivative.

Proposition 2.3 shows that the uniqueness result of Theorem 2.1 is not true without the assumption (2.10) on the spatial derivatives. This suggests that measurement (2.9) alone cannot be used for reconstructing µ. In Section 4.3, working with the same samples µ k as those discussed above, we obtained approximations µ * k of µ k by minimizing a new function H µ , which measures the distance between u(•, x 0 ) and ũ(•, x 0 ). The average relative L 2 -error between µ k and µ k * was 14 times larger than the average relative L 2 -error separating µ k and µ * k . This confirms the usefulness of the spatial derivative measurements for the reconstruction of µ.

Proposition 2 . 3 .

 23 Let µ ∈ M and γ > 0. Assume that α 1 = α 2 and β 1 = β 2 and that u i is symmetric with respect to x = (a+b)/2. Let μ := µ(b-(x-a)) for x ∈ [a, b]. Then, the solutions u and ũ to (P µ,γ ) and (P μ,γ ) satisfy u(t, a+b

. 14 )

 14 Then, we have µ ≡ μ on [a, b] and γ = γ. Consequently u ≡ ũ in [0, +∞) × [a, b]. If β 1 > 0 (resp. β 2 > 0), this statement remains true for x 0 = a (resp. x 0 = b).

  The subset M of C 0,η ([a, b]) made of piecewise analytic functions is much larger than the set of analytic functions on[a, b]. It indeed contains some functions whose regularity is not higher than C 0,η , and some functions which are constant on some subsets of [a, b]. Our results hold true if M is replaced by any subset M ′ of C 0,η ([a, b]) such that for any couple of elements in M ′ , the subset of[a, b] where these two elements intersect has a finite number of connected components.

. 16 )

 16 Proof of Theorem 2.1: In that case γ = γ. Equation (3.15) then reduces to ∂U ∂t -D ∂ 2 U ∂x 2 = μU -γU (u + ũ) + u m.

  iv) b 1 = b, and m = 0 on [x 0 , b]. Assume (i). Then, by definition of b 1 , there exists a decreasing sequence y k → b 1 , y k > b 1 , such that |m(y k )| > 0 for all k ≥ 0. Assume that it exists k 0 such that |m(x)| > 0 for all x ∈ (b 1 , y k 0 ). By continuity, m does not change sign in (b 1 , y k 0 ), and therefore in [x 0 , y k 0 ]. This contradicts the definition of b 1 . Thus, for all k, it exists z k ∈ (b 1 , y k ) such that m(z k ) = 0. (3.18)

  x 0 ] s.t. m has a constant sign on [x, x 0 ]}, the same argument as above shows that a 1 = a and m = 0 on [a, x 0 ]. Thus, finally, m ≡ 0 on [a, b] and this concludes the proof of Theorem 2.1. Proof of Theorem 2.4: From the assumptions (2.12) and (2.14) of Theorem 2.4, equation (3.15) at x = x 0 reduces to

Figure 1 :

 1 Figure 1: (a) An example of function µ in E, together with a function µ * which was obtained by minimizing G µ . In this case µµ * L 2 ([0,1]) / µ L 2 ([0,1]) = 0.03 and G µ (µ * ) = 2 • 10 -6 . (b) The same function µ together with the function µ * obtained by minimizing H µ . Here, µµ * L 2 ([0,1]) / µ L 2 ([0,1]) = 0.47 and H µ (µ * ) = 2 • 10 -6 .

MATLAB's fminunc medium-scale optimization algorithm uses a Quasi-Newton method with a mixed quadratic and cubic line search procedure. Our stopping criterion was based on the maximum number of evaluations of the function Gµ, which was set at

• 10

.

Acknowledgements

The authors would like to thank two anonymous referees for their valuable comments on an earlier version of this paper. The first author is supported by the French "Agence Nationale de la Recherche" within the projects Colon-SGS, PREFERED, and URTICLIM.

Appendix A: maximum principle

The following version of the parabolic maximum principle can be found in [START_REF] Friedman | Partial differential equations of parabolic type[END_REF]Ch. 2] and [START_REF] Protter | Maximum Principles in Differential Equations[END_REF]Ch. 3].

(c) If u ≥ 0, the results (a) and (b) remain true without the assumption c(t, x) ≤ 0.

An immediate corollary of this theorem is:

Proof of Corollary 5.2: Assume that it exists (t * , x * ) ∈ (0, +∞)×(a, b) such that u(t * , x * ) < 0.

Set w(t, x) = u e -λt , for λ > 0 large enough such that

The function w satisfies:

Since w(0, x) = u i (x) ≥ 0 in (a, b) and w(t * , x * ) < 0, the function w admits a minimum m * < 0 in (0, t * ]×[a, b]. From Theorem 5.1 (a), and since u i ≡ 0, this minimum is attained at a boundary point: it exits t ′ ∈ (0, t * ] such that w(t ′ , a) = m * < 0 or w(t ′ , b) = m * < 0. Without loss of generality, we can assume in the sequel that w(t ′ , a) = m * < 0. From Theorem 5.1 (b), we obtain ∂w ∂x (t ′ , a) > 0. Using the boundary conditions in problem (P µ,γ ), we finally get:

Using assumption (2.5), we get a contradiction. Thus u(t, x) ≥ 0 in (0, +∞)× (a, b). The conclusion then follows from Theorem 5.1 (c).

Appendix B: numerical solutions of (P µ,γ ) and (P μ,γ )

The equations (P µ,γ ) and (P μ,γ ) were solved using Comsol Multiphysics time-dependent solver, using second order finite element method (FEM) with 960 elements. This solver uses a method of lines approach incorporating variable order variable stepsize backward differentiation formulas. Nonlinearities are treated using a Newton's method. The interested reader can get more information in Comsol Multiphysics user's guide.


Aurélien Hazan
email: aurelien.hazan@univ-paris1.fr

Michel Verleysen
email: michel.verleysen@uclouvain.be

Marie Cottrell

Jérôme Lacaille
email: jerome.lacaille@snecma.fr

Trajectory Clustering for Vibration Detection in Aircraft Engines

The automatic detection of the vibration signature of rotating parts of an aircraft engine is considered. This paper introduces an algorithm that takes into account the variation over time of the level of detection of orders, i.e. vibrations ate multiples of the rotating speed. The detection level over time at a specific order are gathered in a socalled trajectory. It is shown that clustering the trajectories to classify them into detected and non-detected orders improves the robustness to noise and other external conditions, compared to a traditional statistical signal detection by an hypothesis test. The algorithms are illustrated in real aircraft engine data.

Introduction

We tackle the issue of monitoring the behavior of an aircraft engine from the point of view of measured vibrations. Indeed, abnormal level or odd pattern of vibrations may be the consequence of mechanical or sensor malfunction, both of dramatic importance for engine manufacturers and airline operators.

More precisely this work focuses on the detection of the signature of specific parts of a turbofan engine (the fans), whose vibration levels are modelled as vibration trajectories, i.e. as the amount of vibration of a given fan with respect to time. The trajectories are then clustered thanks to standard clustering algorithms. The obtained clusters gather all trajectories that correspond to fans whose vibratory signature is present in the data.

Section 2 gives a more thorough introduction to the problem. In Section 4 the algorithms are presented. Results are discussed in Section 5, while conclusion and perspectives are sketched in Section 6.

Problem description

A turbofan whose structure is presented by Fig. 1 is considered. Air from the outside enters an intake, then is successively compressed by the low-pressure (LP) and high-pressure (HP) compressors. Compressed air passes to a combustion chamber, where it is mixed with fuel and burnt. Both compressors are powered by turbines located at the rear of the engine, which transmit their energy to the compressors through two contra-rotating shafts, the low-pressure (LP) shaft and the high-pressure (HP) shaft.

Although turbofan condition monitoring can be achieved in various ways, we take the stand to focus on vibrations monitoring in this work. Two accelerometers provide vibration measurements at a constant 51 kHz frequency. Since compressors and turbines are fan-like components made of a varying number of blades mounted on the shafts, it is expected that their motion entails vibrations at frequencies which are multiples of shaft speeds. Although this is a strong simplification of the overall vibratory behavior of a bladed disk mounted on a rotor [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling[END_REF], it allows an efficient detection of malfunction and damages, given the low quantity of available information (two vibration sensors). Vibration patterns corresponding to multiples of shaft speed are known as "orders" in the engineering field. Moreover, vibrations signals are usually processed not in the time-domain, but in the frequency [START_REF] Braun | Mechanical Signature Analysis: theory and Applications[END_REF] or in the time-frequency domain [3, 15.6]. When the rotation speed of the engine is constant, vibration signals can be considered as stationary, and the classical Fourier transform is sufficient. However when the rotation speed is varying, signals are non-stationary and the spectrogram or more advanced time-frequency distributions might be of some help. Examples of a Fourier transform and a spectrogram on a vibration signal are given in Fig. 2(a,b). In addition to accelerometers, the sensors give an approximation of both shaft velocities. Since many mechanical parts of the engine rotate at a speed multiple of these shafts speeds (compressors and turbines for example), we can anticipate the position of their signature in the frequency or time-frequency domain, as illustrated in Fig. 2(c) where an unknown signature in the spectrogram and an expected signature superpose at a particular order.

As we can notice from Fig. 2(b), the signature is embedded in noise. For specific orders, it remains clearly visible whereas for others it vanishes. Furthermore, there is a potentially large number of orders. For both reasons, order detection is most of the time the task of an expert. Automatic detection would hence accelerate the process, and relieve the expert of this daunting task.

Several difficulties arise when considering the issue of automatic detection of orders in vibration data. First, as already mentioned, the signal is noisy. Second, the magnitude of vibration of the orders is a function of the excitation. The latter is a direct function of the shafts speeds. Indeed, the fans behave as mechanical parts that are excited and resonate for specific frequencies only. Consequently, if we examine a sufficiently long recording, the orders will appear with varying magnitude along time. The frequency response of rotating machines is classically summarized in Campbell diagrams where the experimentally4 measured vibration response spectrum appear as a function of the rotation speed of the shaft. If we suppose that the underlying shaft rotation speed follows a regular enough time pattern, then there should be some degree of continuity in the vibration magnitude of the orders. Admitting that signals are split into a finite number of windows for computational reasons, we propose to consider vibration magnitude trajectories, where to each time window corresponds a robust measure of the vibration magnitude at a given order, for a selected shaft. This needs to be done systematically since we ignore whether a given order k of LP or HP shaft is visible for a specific regime. At the end, a large number of trajectories will be available, which we propose to cluster thanks to dedicated algorithms in order to segregate those who correspond to orders that are present in the data from those who do not. This is summarized by Fig. 3. The expected outcome of this algorithm is thus a cluster of trajectories that correspond to what the expert would have singled-out as orders that are actually present in the recordings.

Related works are given a quick review in Section 3, while a detailed presentation of available data is given in Section 4.1. The proposed method is discussed in Section 4.2. Canonical signal processing detection procedures in the case of partially unknown sine wave are evoked in Section 4.3 for benchmarking purposes.

Related work

This work addresses the problem of feature detection, where the features are the vibration signatures of bladed disks that compose compressors and turbines of an aircraft engine.

It is intended to be part of a Condition Based Monitoring (CBM) framework. CBM for industrial machines has been attracting increasing attention over the years in both academic and industrial areas. According to [START_REF] Randall | State of art in monitoring rotating machinery -part I[END_REF] it consists in four main steps: data acquisition, feature extraction, feature selection, and decisionmaking. The first two steps rely on mechanical modeling or rotor dynamics [START_REF] Muszynska | Rotordynamics. Taylor and Francis[END_REF], noise and vibration phenomena in rotating machines [START_REF] Lyon | Machinery Noise and Diagnostics[END_REF][START_REF] Randall | State of art in monitoring rotating machinery -part I[END_REF] and data analysis [START_REF] Peng | Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography[END_REF]. The latter builds on the general tools and methods developped in signal processing [START_REF] Mallat | Une exploration des signaux en ondelettes[END_REF][START_REF] Boashash | Time-frequency signal analysis and processing -A comprehensive reference[END_REF], statistical signal estimation and detection [START_REF] Kay | Fundamentals of statistical signal processing: detection theory[END_REF][START_REF] Poor | An introduction to signal detection and estimation 2nd ed[END_REF][START_REF] Van Trees | Detection, estimation, and modulation theory Part 1[END_REF], learning theory, change detection [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF], fault detection and isolation [START_REF] Gertler | Fault detection and diagnosis in engineering systems[END_REF].

Aircraft CBM deals with many problems such as structural health monitoring. It treats engine health monitoring (EHM) as a special case [START_REF] Tumer | A survey of aircraft engine health monitoring systems[END_REF]15]. As a subtopic of EHM, vibrations monitoring in engines addresses the following issues: rotor/stator contact [START_REF] Peng | Detection of the rubbing-caused impacts for rotor-stator fault diagnosis using reassigned scalogram[END_REF], rotor unbalance, blade defects [START_REF] Kharyton | fault detection of blades in blades ov an aviation engines in operation[END_REF], bearing [START_REF] Orsagh | Prognostics/diagnostics for gas turbine engine bearings[END_REF] and gearings defects [19]. Another important topic is order tracking, i.e. the precise estimation of the frequency and amplitude of a periodic signal whose leading frequency is varying in time. This may involve data resampling in the case of rotating machines [START_REF] Potter | Computed order tracking obsoletes older methods[END_REF][START_REF] Fyfe | Analysis of computed order tracking[END_REF], non-parametric time frequency methods such as the Gabor transform [START_REF] Qian | Gabor expansion for order tracking[END_REF] or parametric methods such as the Vold-Kalman filter [START_REF] Vold | Resolution order tracking at extreme slow rates, using Kalman tracking filters[END_REF]24], that models the vibrating machine and use estimation tools to track the frequencies and amplitudes of orders.

When features have been extracted, decisions can be made, such as change (or novelty) detection. The decision concerning the health of the engine is either taken from a statistical viewpoint [START_REF] Basseville | Analyse et surveillance vibratoire d'une machine en rotation[END_REF][START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF], or can derive from learning techniques such as neural networks [START_REF] Ypma | Learning methods for machine vibration analysis and health monitoring[END_REF][START_REF] Staszewski | Signal processing for damage detection[END_REF].

In this work the aim is not to detect an anomaly, but a normal bladed disk signal. Up to our knowledge, there is no example of such order detection in the litterature. The case of novelty or abrupt change detection in engine have been covered, as well as feature extraction for many particular mechanical parts, but not for bladed disks in turbines and compressors. The works dedicated to these parts aim at estimating their rotation speed, amplitude and phase, but not to detect their presence. In the following we show that classical tools such as sine wave detection and unsupervised clustering enable us to tackle this issue.

Methods and data

Data

The recordings under study were provided by the Health Monitoring Department of SNECMA 5 and correspond to a dual-shaft turbofan mounted on a testbench, that undergoes a continuous acceleration during several minutes. They include raw vibration outputs of an accelerometer sampled at 51kHz, as well as LP and HP shaft angular velocity computed from raw keyphasor data and sampled at 6.25 Hz. Sample time series are plotted in Fig. 4.

Clustering contrast trajectories for order detection

First we give a precise meaning to what was termed "vibration magnitude trajectory" in Section 2. Conforming to current practice in the field of vibration monitoring for rotating machine, we choose a classical time-frequency representation, the Discrete Gabor Transform (DGT) [START_REF] Qian | Gabor expansion for order tracking[END_REF], which is a special kind of Discrete Short-Time Fourier Transform [START_REF] Mallat | Une exploration des signaux en ondelettes[END_REF][START_REF] Feichtinger | Gabor analysis and algorithms: theory and applications[END_REF]. The representation of a given vibratory signal is a matrix, indexed in both time and frequency.

As stated earlier, we focus on the signature of "orders", produced by the bladed disks that compose the compressors and turbines. Their simplified signature in the time-frequency domain has been illustrated by Fig. 2(c). We write G(t, ν) the complex-valued time-frequency transform at discrete time t and discrete frequency ν. Since we know approximately the frequency of an order from the choice of an integer number (the order) and the estimation of the shaft speed, we can build a region in the time-frequency plane in which we expect the order signature to appear, if it is present in vibratory data. This operation is called masking and is illustrated in Fig. 5(a).

Then, having focused on a specific order, we still need to compute a measure of magnitude for that order. Several statistical preprocessing steps (computation of the modulus, of the logarithm, then centering and scaling of the data) are followed by a geometric flattening of the signature as shown in Fig. 5(b). We call G ′ (t, ν) the resulting quantity. Finally the flattened signature is time-average, then the contrast function is computed.

∀ν, G(ν) = 1 T t2 t=t1 G ′ (t, ν) (1)
c k (n) = max ν G(ν) -min ν G(ν) (2)
where k and n are respectively the order and the window indices, T is the length of the time window, and t 1 and t 2 are respectively the start and end of that time window. Secondly, we describe the trajectory clustering algorithm announced in 2. The previous masking step returns a contrast value for each order k in a userdefined range [k 1 , k 2 ], and for each window n ≤ N that divides the full vibration recording. Hence for each order k we get an ordered set of scalar real-valued contrasts T k = {c k (1), . . . , c k (N )} indexed by the order k. This set T is what we call a trajectory. We then propose to cluster it into a fixed number C = 3 of clusters, with the k-means algorithm. The whole procedure is illustrated in Fig. 6. Clustering trajectories is a way to classify them into classes, corresponding respectively to orders that are detected, not detected, and a possible intermediate, ambiguous classes. One could argue that if an order k is effectively present in the recordings, any contrast c k (n) belonging to T k should be high, allowing a detection on each time window n. However, one has to take into account the fact that the vibratory signals are extremely noisy. Splitting into windows is thus a way to expect an improved detection robustness in some windows. Averaging the detection hits over the windows could then improve the detection results, compared to a single detection (over one window or the whole signal).

Still, this average would completely lose the temporal correlation between orders. Indeed it is naturally expected that at some engine speeds (corresponding to some of the windows), all or most of the orders will be hardly detected (either because a large part of the engine will be less subject to vibrations, or because the latter will be polluted by higher noise). It is therefore of primary importance to take into account the longitudinal aspect of the detections at various orders, i.e. the temporal correlation. This is the justification for the original use of trajectories in this work.

In Section 4.3 we summarize classical results relevant to our problem in the field of statistical signal detection, in order to assess the performance of the proposed method.

Statistical signal detection method

From the point of view of classical statistical signal detection theory [START_REF] Kay | Fundamentals of statistical signal processing: detection theory[END_REF][START_REF] Poor | An introduction to signal detection and estimation 2nd ed[END_REF][START_REF] Van Trees | Detection, estimation, and modulation theory Part 1[END_REF], the problem formulated in Section 2 is an instance of the detection of a deterministic sine wave in noise, with unknown parameters [9, 7.6.2]. Indeed, with the simplifying hypotheses stated in Section 2, the compressor and turbine disk blades whose signature is under study are seen as curves in the time-frequency plane, i.e. as sine waves with slowly varying leading frequency. On time windows where the frequency can be considered as constant, this sine wave can be parametrized as follows:

s(i) = A cos(2πf i + φ) (3)
where A is the amplitude, i the time instant, f the frequency and φ the phase difference. Apart from this wave the signal contains a vast amount of other deterministic or random contributions which we call "noise" since they are irrelevant to the purpose of detecting a given order of a specific shaft. The detection of the sine wave in the noise is formalized as the following hypothesis test:

H 0 : z20(i) = w(i) ∀i ∈ [0, I -1] H 1 : z20(i) = w(i) + A cos(2πf i + φ) ∀i ∈ [0, I -1]
where z20 is the accelerometric signal, w(i) is the noise term. It should be noted that, while the frequency f is approximately known, the amplitude A and phase difference φ are unknown to the observer. Under noise normality and independence hypotheses that diserve further discussion, this hypothesis test is amenable to analytic treatment under the Generalized Likelihood Ratio Test framework. Since this material is standard, we merely state the corresponding results, which are given proper development in classical textbooks (see for example [9, 7.6

.2]) :

the decision rule is written I(f 0 ) ≷ H1 H0 γ ′ , where the decision statistics is I(f 0 ), the value of the periodogram evaluated at target frequency f 0 which can be approximated by Discrete Fourier Transform. γ ′ is the detection threshold.

the false alarm probability is α = expγ ′ σ 2 , where σ is the standard deviation of the noise, supposedly known. Usually α takes a fixed user-defined value such as α = 1% so that the value of the detection threshold γ ′ can be easily derived.

the performance of this test is quantified by the power π of the test, whose analytic expression can be computed.

With these results, we can build a sine wave detector. Given a vibratory time series and a few contextual parameters such as the noise variance, the detector will accept or reject hypothesis H 0 depending if the signal is composed of noise or if an order is present in the signal. Note that several improvements over this standard procedure should be made, because the estimated noise does not respect exactly the white noise hypothesis. For example, the energy is not homogenous in the frequency spectrum. Consequently we divide the frequency spectrum into bins and perform the above detection separately in each bin. In Section 5 we compare the results with those elicited by the clustering algorithm.

Results

We first illustrate the computation of the contrast function defined in Section 4.2.

The first steps depicted by Fig. 5 are the statistical and geometric preprocessing steps. Fig. 7(a) shows the resulting preprocessed time-frequency transforms for six chosen orders of the HP shaft. Three orders (38,53,68) correspond to the number of blades of three bladed disks belonging to the compressor, therefore we expect the corresponding signatures to be present in vibration signals. Three ordrers (20,90,100) were selected because no significant vibratory activity is expected. In Fig. 7(a) thick lines appear clearly in the first three cases, whereas no specific pattern except the background noise can be noticed in the last three cases. Time averaged values of the time-frequency transforms are plotted in Fig. 7(b), as well as their peak value in Fig. 7(c). The latter clearly shows that peak values are higher for orders 38, 53, 68 that correspond to actual mechanical rotating parts.

Secondly we comment on the contrast trajectory clusters that are found out by the clustering algorithm mentioned in Fig. 6. The number of clusters was set to C = 3. Results are shown in Fig. 8(a). The first cluster gathers trajectories with high mean value over the window range, that increases at the beginning, decreases at the end, and experiences a sudden decrease at window index 6. In the second one, trajectories are flat and have low mean value. The third cluster mixes different types of trajectories that have an average intermediate value.

In first approximation, only the trajectories in cluster 1 are meaningful to the detection task.

In order to assess the significance of results, we use prior mechanical information. Indeed, the number of blades that compose the compressor and turbine mounted on HP shaft is known. Because of many mechanical factors it is not certain whether or not each bladed disk will have a noticeable vibratory activity, but this information can be helpful for comparison purposes. In Fig. 8(b), we plot the composition of the clusters obtained above (indices 1 and 3), and compare it with the cluster built with prior information which is given cluster label 0. Lastly, clusters found by the statistical detection algorithm from Section 4.3 are labelled as cluster D. We see that between orders 35 to 80 many orders are shared by clusters 0 and 1. This is true also for clusters 0 and D. However, clusters 2 and 3 show little similarity with cluster 0. This is coherent with our initial expectation, stated at the end of Section 2. In addition we remark that many low orders are detected in cluster D. This could be explained by the fact that low frequencies bear more energy that higher ones, as evidenced in Fig. 2(b). The clustering method is less prone to overweigthing such orders because of the higher energy content in the low orders area. Lastly the fact that orders not expected from mechanical knowledge appear both in clusters 1 and D suggest that interesting information not provided by naive mechanical data was actually discovered.

Lastly we confirm quantitatively these observations. Mutual information [START_REF] Kraskov | Estimating mutual information[END_REF] is used here to measure the similarity between clusters. Setting cluster 0 as the default one, we compute the mutual information between the corresponding column of the composition matrix displayed in Fig. 8(b) and the other columns. What we expect is that self-information (i.e. the information between column 0 and itself) is high, while information between unrelated clusters such as 0 and 2 should be low. Moreover, mutual information between cluster pairs 0/1 and 0/D should be between the self-information value, and the unrelated clusters value. Indeed this is what we observe from Fig. 9, where the highest value corresponds to cluster pairs 0/0, while cluster pair 0/1 and 0/D have high mutual information. The negative value, which is theoretically meaningless, is a known limitation of the estimation algorithm. Lastly, the mutual information for cluster pair 0/D is lower than for cluster pair 0/1, which seems to indicate a higher performance of the proposed algoritm. The good results obtained by unsupervised clustering over signal-processing algorithms are suprising at first mainly because additional knowledge is embedded in the latter. Indeed, one needs to model the signal and the noise components before deriving a decision algorithm. We suggest that the results of the signal detection algorithm could be explained first by poor agreement between real data and noise hypotheses. This should be examined in further experiments.

Nevertheless, it remains that in the clustering approach we take into account the continuity of contrast values in time, which is a consequence of the continuity of shafts rotation speeds as a function of time. This continuity remains unexploited by the classical signal-detection algorithm, which iterates the decision process over successive windows without relating them.

Conclusion and perspectives

In this work we have tackled the issue of order detection, i.e. the discovery of vibration patterns in noisy aircraft engine vibration signals. We proposed a contrast measure whose aim is to single out significant vibration patterns that correspond to compressor and turbines mounted on the engine shafts. Then a clustering algorithm is built, and compared to a statistical signal detection procedure. We show with real data that the clustering method performs well, and give quantitative measure of this performance. Future works will aim at :

increasing the statistical significance by enlarging the database to several engines, in both acceleration and deceleration situations. Theoretical properties of the estimators of contrast, and of mutual information could be studied as well.

improving the signal detection algorithm, for example by considering extension to colored noise situation. In addition the continuity from one window to the following could be used. -merging the decisions from both methods.

refining the clustering, mainly the interpretation of the intermediate cluster.

assessing the continuity hypothesis and using it as prior for clustering, from the knowledge of theoretical Campbell diagrams.

Fig. 1 .

 1 Fig. 1. Turbofan engine. Simplified diagram of fan, low-pressure and high-pressure compressors and turbines attached to their respective shafts.

Fig. 2 .

 2 Fig. 2. Vibration analysis of accelerometric data : (a) in the frequency domain with the Discrete Fourier Transform; (b) in the time-frequency domain with the spectrogram; (c) idealized signature of fan-like mechanical part in accelerometer data. An unknown signature in the spectogram and an expected signature at 36 times the frequency of LP shaft superpose .

Fig. 3 .

 3 Fig. 3. Idealized diagram of trajectories clustering.

Fig. 4 .

 4 Fig. 4. (a) Raw accelerometric data; (b) LP shaft angular velocity in rpm.

Fig. 5 .

 5 Fig. 5. (a) Order masking takes as its inputs the vibratory signal, a specific order integer k and the shaft rotation speed N1 or N2. The expected main frequency is f k,N 1 , and the region that surrounds ths central frequency has a fixed user-defined width ; (b) contrast.

Fig. 6 .

 6 Fig. 6. Clustering algorithm.

Fig. 7 .

 7 Fig. 7. Contrasts of orders 38,53,68,20,90,100 of HP-shaft N2: (a) DGT after preprocessing; (b) mean DGT along time axis; (c) contrasts.

Fig. 8 .

 8 Fig. 8. Comparison of clustering and classical detection results for z20 vibratory data: (a) contrast vectors clustered by k-means with C = 3 clusters; (b) composition matrix for each cluster. Dark line with matrix coordinate i, j means that order i belongs to cluster j. Column 0 represents orders for which fan signature is expected from mechanical knowledge. Columns 1 to 3 are produced by the clustering algorithm. Column D is output by the statistical signal detection algorithm.

Fig. 9 .

 9 Fig. 9. Mutual information between pairs of clusters.

such diagrams can be drawn from analytical models as well, if the eigenfrequencies can be computed.

http://www.snecma.fr

Acknowledgements

We thank division YY of Snecma for providing us with vibrations data, more particularly S. Blanchard and J. Griffaton. We also thank E. Côme for suggesting to compute similarity measure between clusters with mutual information. Discrete Gabor Transform is computed thanks to the LTFAT library [30], available at http://ltfat.sourceforge.net/. Mutual information is computed thanks to code made publicly available by Kraskov at http://www.klab.caltech.edu/ kraskov/MILCA/.