D Y Liu

T M Laleg-Kirati

O Gibaru

W Perruquetti
email: wilfrid.perruquetti@inria.fr

Identification of fractional order systems using modulating functions method

Da-Yan Liu, Taous-Meriem Laleg-Kirati, Olivier Gibaru and Wilfrid Perruquetti

Abstract-The modulating functions method has been used for the identification of linear and nonlinear systems. In this paper, we generalize this method to the on-line identification of fractional order systems based on the Riemann-Liouville fractional derivatives. First, a new fractional integration by parts formula involving the fractional derivative of a modulating function is given. Then, we apply this formula to a fractional order system, for which the fractional derivatives of the input and the output can be transferred into the ones of the modulating functions. By choosing a set of modulating functions, a linear system of algebraic equations is obtained. Hence, the unknown parameters of a fractional order system can be estimated by solving a linear system. Using this method, we do not need any initial values which are usually unknown and not equal to zero. Also we do not need to estimate the fractional derivatives of noisy output. Moreover, it is shown that the proposed estimators are robust against high frequency sinusoidal noises and the ones due to a class of stochastic processes. Finally, the efficiency and the stability of the proposed method is confirmed by some numerical simulations.

I. INTRODUCTION

Fractional differential equations and fractional integrals are gaining importance in research community because of their capacity to accurately describe real world processes. The flow of fluid in a porous media, the conduction of heat in a semiinfinite slab, the voltage-current relation in a semi infinite transmission line are such examples of processes naturally modeled by fractional differential equations or fractional integrals.

This paper is dealing with the identification of fractional order dynamical systems. The identification of such systems has been used for instance, for the estimation of the state of charge of lead acid batteries [START_REF] Sabatier | Fractional system identification for lead acid battery state of charge estimation[END_REF], and for the identification in thermal systems [START_REF] Gabano | Fractional modeling and identification of thermal systems[END_REF], [START_REF] Gabano | Identification of a thermal system using continuous linear parameter-varying fractional modelling[END_REF]. The goal of system identification is to estimate the parameters of a model from system input/output measurements. Different methods have been proposed for the identification of fractional order systems. Most of them consist in the generalization to fractional order systems of standard methods that were used in the identification of systems with integer order derivatives. We can classify these methods into time domain methods and frequency domain methods. Timedomain methods have been introduced for example in [START_REF] Oustaloup | Identification of non integer order system in the time domain[END_REF], [START_REF] Djouambia | Recursive prediction error identification of fractional order models[END_REF], where a method based on the discretization of a fractional differential equation using Grunwald definition has been introduced, and the parameters have been estimated using least square approach. In [START_REF] Trigeassou | Modeling and identification of a non integer order system[END_REF], a method based on the approximation of a fractional integrator by a rational model has been proposed. In [START_REF] Aoun | Synthesis of fractional Laguerre basis for system approximation[END_REF], the use of methods based on fractional orthogonal bases has been introduced. Other techniques can be also found for example in [START_REF] Cois | Fractional state variable filter for system identification by fractional model[END_REF], [START_REF] Cois | Non integer model from modal decomposition for time domain system identification[END_REF], [START_REF] Malti | Tutorial on system identification using fractional differentiation models[END_REF] and the references therein.

In this paper, we are interested in the identification of fractional order systems using modulating functions method in case of noisy measurements. Modulating functions method has been developed by Shinbrot in [START_REF] Shinbrot | On the analysis of linear and nonlinear dynamic systems from transient-response data[END_REF], [START_REF] Shinbrot | On the analysis of linear and nonlinear systems[END_REF] to estimate the parameters of a state space representation. Thanks to the properties of modulating functions, the fractional differential equation defining a fractional order system is transformed into a linear system of algebraic equations. Hence, instead of solving a fractional differential equation where the initial values are often unknown, the problem of identification is transformed into solving a linear system where the initial conditions are not required. Generalization of modulating functions method to fractional systems has been already proposed, for example in [START_REF] Janiczek | Generalization of modulating functions method in the fractional differential equations[END_REF]. However, the authors of this paper proposed only to reduce the orders of the derivatives in a fractional differential equation. Moreover, the noisy case has not been considered.

In the next section, basic definitions of fractional derivatives and modulating functions are recalled. Then in Section III, modulating functions method is applied to the identification of fractional order linear systems. Error analysis in both the continuous and discrete cases are presented in Section IV. Numerical results are presented in Section V, followed by conclusions summarizing the main results obtained.

II. PRELIMINARY

In this section, first we recall the definitions and some useful properties of the Riemann-Liouville fractional derivative and modulating functions. Then, a new fractional integration by parts formula is given.

A. Riemann-Liouville fractional derivative

Let f be a continuous function defined on R, then the Riemann-Liouville fractional derivative of f is defined by (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] p. 62):

∀ t ∈ R * + , D α t f (t) := 1 Γ(l -α) d l dt l t 0 (t -τ ) l-α-1 f (τ ) dτ, (1)
where l -1 ≤ α < l with l ∈ N * , and Γ(•) is the Gamma function (see [START_REF] Abramowitz | editeurs. Handbook of mathematical functions[END_REF] p. 255). As an example, using (1) the fractional derivative of an n th (n ∈ N) order polynomial is given by (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] p. 72):

∀ t ∈ R * + , D α t t n = Γ(n + 1) Γ(n + 1 -α) t n-α . (2)
We assume that the fractional derivative and the Laplace transform of f both exist, then the Laplace transform of the fractional derivative of f is given by (see [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF] p. 284): ∀ s ∈ C,

L {D α t f (t)} (s) = s α f (s) - l-1 i=0 s i D α-i-1 t f (t) t=0 , (3)
where f denotes the Laplace transform of f , and s denotes the variable in the frequency domain.

Finally, we recall some results on the existence and the initial values of the fractional derivative in the following proposition.

Proposition 1: (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] pp. 75-77

) If f ∈ C l-1 ([0, T ]) and f (l) ∈ L([0, T ]) where T ∈ R * + , then the Riemann-Liouville fractional derivative D α t f (t) exists, where l -1 ≤ α < l with l ∈ N * . Moreover, the condition f (i) (0) = 0, for i = 0, . . . , l- 1, is equivalent to the following condition: [D α t f (t)] t=0 = 0.

B. Modulating functions

Let l ∈ N * , T ∈ R * + , and g be a function satisfying the following properties:

(P 1 ) : g ∈ C l ([0, T ]); (P 2 ) : g (j) (0) = g (j) (T ) = 0, ∀ j = 0, 1, .
. . , l -1, then g is called modulating function of order l [START_REF] Preising | Theory and application of the modulating function method. I: Review and theory of the method and theory of the splinetype modulating functions[END_REF]. Hence, according to Proposition 1, the modulating function g has also the following properties:

(P 3 ) : ∀ 0 ≤ β < l, D β t g(t) exists; (P 4 ) : ∀ 0 ≤ β < l, D β t g(t) t=0 = 0.

C. Fractional integration by parts

Another useful property of the modulating functions is given in the following theorem.

Theorem 1: Let y be a function such that the α th order fractional derivative exits and g be a modulating function of order l with l -1 ≤ α < l with l ∈ N * . Then, we have:

T 0 g(T -t) D α t y(t) dt = T 0 D α t g(t) y(T -t) dt, (4)
where T ∈ R * + . Proof: By applying the convolution theorem of the Laplace transform (see [START_REF] Abramowitz | editeurs. Handbook of mathematical functions[END_REF], p. 1020), we get:

L T 0 g(T -t) D α t y(t) dt (s) = ĝ(s)L {D α t y(t)} (s).
(5) Then, using (3) we obtain:

ĝ(s)L {D α t y(t)} (s) = ĝ(s)s α ŷ(s) - l-1 i=0 s i ĝ(s) D α-i-1 t y(t) t=0 . (6)
Moreover, using ( 3) and (P 4 ) we get:

L {D α t g(t)} (s) = s α ĝ(s), (7)
L g (i) (t) (s) = s i ĝ(s), (8)
for i = 0, . . . , l -1. Consequently, by applying ( 7), ( 8) and the inverse of the Laplace transform to (6), we obtain:

L -1 {ĝ(s)L {D α t y(t)} (s)} (T ) = L -1 {L {D α t g(t)} (s) ŷ(s)} (T ) - l-1 i=0 g (i) (T ) D α-i-1 t y(t) t=0 . (9)
Using (P 2 ), the initial conditions D α-i-1 t y(t) t=0 , for i = 0, 1, . . . , l -1, can be eliminated.

Finally, this proof can be completed by applying the convolution theorem of the Laplace transform to [START_REF] Cois | Non integer model from modal decomposition for time domain system identification[END_REF].

According to the previous theorem, we can see that by working in the frequency domain: on the one hand, we can obtain an integral formula which can be considered as the generalization of the classical integration by parts formula. Another fractional integration by parts formula has also been given in [START_REF] Podlubny | Adjoint fractional differential expressions and operators[END_REF], however the Caputo fractional derivative was involved in the formula; on the other hand, the initial conditions of the fractional derivatives of y can be eliminated using a modulating function. In fact, the idea of obtaining this theorem is inspired by the recent algebraic parametric estimation technique [START_REF] Fliess | An algebraic framework for linear identification[END_REF], [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF], [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF], [START_REF] Liu | Parameters estimation of a noisy sinusoidal signal with time-varying amplitude[END_REF], [START_REF] Liu | An error analysis in the algebraic estimation of a noisy sinusoidal signal[END_REF], [START_REF] Ushirobira | Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data[END_REF], [START_REF] Perruquetti | An algebraic approach for Humane posture estimation in the sagital plane using accelerometer noisy signal[END_REF], [START_REF] Trapero | An algebraic frequency estimator for a biased and noisy sinusoidal signal[END_REF], which eliminates the unknown initial conditions by applying algebraic manipulations in the frequency domain.

III. IDENTIFICATION OF FRACTIONAL ORDER SYSTEMS

A. Fractional order linear systems

In this section, we consider a class of fractional order linear systems which are defined by the following fractional differential equation:

∀ t ∈ I = [0, T ] ⊂ R * + , L i=0 a i D αi t y(t) = M j=0 b j D βj t u(t), ( 10
)
where y is the output, u is the input, a i , b j ∈ R * are unknown parameters to be identified, and

α i , β j ∈ R + are assumed 0 ≤ α 0 < α 1 < • • • < α L , 0 ≤ β 0 < β 1 < • • • < β M with L, M ∈ N.
Let y ̟ be a noisy observation of y on the interval I:

∀ t ∈ I, y ̟ (t) = y(t) + ̟(t), (11)
where ̟ is an integrable noise 1 . We are going to estimate the unknown parameters in (10) using the input u and the observation y ̟ . One of the standard methods for the identification of systems with integer order derivatives is to use the leastsquares method [START_REF] Ljung | System Identification: Theory for the User[END_REF]. This method was generalized to fractional order systems in [START_REF] Oustaloup | Identification of non integer order system in the time domain[END_REF], [START_REF] Djouambia | Recursive prediction error identification of fractional order models[END_REF]. In this method, we need to estimate the fractional derivatives of y and u using a fractional order differentiator [START_REF] Chen | Digital Fractional Order Savitzky-Golay Differentiator[END_REF], [START_REF] Liu | Fractional order differentiation by integration with Jacobi polynomials[END_REF], [START_REF] Liu | Non-asymptotic fractional order differentiators via an algebraic parametric method[END_REF]. However, a fractional order differentiator often contains a truncated term error. An alternative method consists in solving the problem in the frequency domain by applying the Laplace transform to [START_REF] Malti | Tutorial on system identification using fractional differentiation models[END_REF]. However, according to (3), this application can produce unknown initial conditions. In the next subsection, we are going to apply modulating functions method to eliminate these unknown initial conditions.

B. Application of modulating functions method

We denote W = L + M + 1, α = max (α L , β M ), and l = ⌈α⌉, where ⌈α⌉ denotes the smallest integer greater than or equal to α. Then, we take a set of modulating functions

{g n } N n=1 with W ≤ N ∈ N.
Using Theorem 1, we can give the following proposition.

Proposition 2: Let {g n } N n=1 be a set of modulating functions of order l. If we assume that b 0 = 1 in the fractional order linear system defined by [START_REF] Malti | Tutorial on system identification using fractional differentiation models[END_REF], then the unknown parameters in this fractional order linear system can be estimated by solving the following linear system:

U N Y ̟ N B Ã = I N , (12)
where

B = b1 , • • • , bM T , Ã = (ã 0 , • • • , ãL )
T are estimators of the unknown parameters, and

U N (n, j) = - T 0 D βj t g n (t) u(T -t) dt, (13)
Y ̟ N (n, i + 1) = T 0 D αi t g n (t) y ̟ (T -t) dt, (14)
I N (n) = T 0 D β0 t g n (t) u(T -t) dt, (15)
for n = 1, . . . , N , j = 1, . . . , M and i = 0, . . . , L. Proof: By multiplying the modulating functions g n to the equation [START_REF] Malti | Tutorial on system identification using fractional differentiation models[END_REF] and by integrating between 0 and T , we get:

T 0 L i=0 a i g n (T -t) D αi t y(t) dt = T 0 M j=0 b j g n (T -t) D βj t u(t) dt, (16)
for n = 1, . . . , N . Then, using Theorem 1, we get

M j=1 b j U N (n, j)+ L i=0 a i Y N (n, i+1) = T 0 D β0 t g n (t)u(T -t)dt, (17) where
U N (n, j) = - T 0 D βj t g n (t) u(T -t) dt, and Y N (n, i + 1) = T 0 D αi t g n (t) y(T -t) dt.
Finally, this proof can be completed by substituting y by y ̟ in [START_REF] Preising | Theory and application of the modulating function method. I: Review and theory of the method and theory of the splinetype modulating functions[END_REF].

Consequently, thanks to Theorem 1, instead of estimating (resp. calculating) the fractional derivatives of y (resp. u), we calculate the fractional derivatives of the modulating functions. On the one hand, comparing to y, the modulating functions are known and without noise. On the other hand, if the fractional derivatives of u cannot be analytically calculated or are difficult to calculate, we can solve the problem by calculating the ones of the modulating functions.

Finally, let us mention that since the proposed estimators are given in causal case, if we take the value of T to be equal to the time where we estimate the parameters, then these estimators can be used for on-line identification applications.

IV. ERROR ANALYSIS

In this section, we are going to study the noise effect in the integrals obtained in Proposition 2. For this purpose, we study the noise error contributions due to a high frequency sinusoidal noise and the ones due to a class of stochastic processes in continuous case and in discrete case, respectively.

A. Error analysis in continuous case

There are many applications where the output signal is corrupted by a sinusoidal noise of higher frequency [START_REF] Trapero | An algebraic frequency estimator for a biased and noisy sinusoidal signal[END_REF]. Hence, we assume that the noise ̟ is a high frequency sinusoidal noise in this subsection.

By writing y ̟ = y + ̟, the integral Y ̟ N (n, i + 1) given in ( 14), for n = 1, . . . , N and i = 0, . . . , L, can be divided into:

Y ̟ N (n, i + 1) = Y N (n, i + 1) + e ̟ N (n, i + 1), (18)
where Y N (n, i + 1) is given in [START_REF] Preising | Theory and application of the modulating function method. I: Review and theory of the method and theory of the splinetype modulating functions[END_REF], and the associated noise error contribution is given by:

e ̟ N (n, i + 1) = T 0 D αi t g n (t) ̟(T -t) dt. (19)
Consequently, the estimation errors for the estimators given in Proposition 2 only come from these noise error contributions.

In the following proposition, we give error bounds for these noise error contributions. Proposition 3: We assume that ∀ t ∈ I, ̟(t) = c sin(ωt + φ) with c, ω ∈ R * + and φ ∈ [0, 2π[. Moreover, we assume that D αi+1 t g n (t) exists and is continuous on [0, T ], for n = 1, . . . , N and i = 0, . . . , L. Then, we have:

|e ̟ N (n, i + 1)| ≤ c ω T C αi+1 + c ω |[D αi t g n (t)] t=T | , (20)
where

C αi+1 = sup t∈[0,T ] D αi+1 t g n (t)
, and e ̟ N (n, i+1) is given by [START_REF] Fliess | An algebraic framework for linear identification[END_REF].

Proof:

Since d dt {D αi t g n (t)} = D αi+1 t
g n (t) exists, then by applying integration by parts and (P 4 ), we get:

e ̟ N (n, i + 1) = T 0 D αi t g n (t) c sin(ω(T -t) + φ) dt = - c ω T 0 D αi+1 t g n (t) cos(ω(T -t) + φ) dt + c ω cos(φ) [D αi t g n (t)] t=T . (21)
If D αi+1 t g n (t) ∈ C([0, T ]), then this proof can be completed using [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF].

According to the previous proposition, if the frequency of the sinusoidal noise is high, then the associated noise error contributions can be negligible. Consequently, the estimators given in Proposition 2 can cope with this kind of noises.

B. Error analysis in discrete case

From now on, we assume that the noisy observation y ̟ defined in [START_REF] Shinbrot | On the analysis of linear and nonlinear dynamic systems from transient-response data[END_REF] is given in a discrete case. Let y ̟ (t j ) = y(t j ) + ̟(t j ) be a noisy discrete observation of y given with an equidistant sampling period T s , where T s = T m , m ∈ N * , and t j = jT s , for j = 0, • • • , m.

Since y ̟ is a discrete measurement, we apply a numerical integration method to approximate the integrals in [START_REF] Shinbrot | On the analysis of linear and nonlinear systems[END_REF]. Let w 0 ≥ 0, w m ≥ 0 and w j > 0 for j = 1, . . . , m -1 be the weights for a given numerical integration method, where weight a 0 (resp. a m ) is set to zero when there is an infinite value at t 0 = 0 (resp. t m = T ). Then, the integral Y ̟ N (n, i+1) given in [START_REF] Podlubny | Fractional Differential Equations[END_REF], for n = 1, . . . , N and i = 0, . . . , L, can be approximated by:

Y ̟ N,m (n, i + 1) := T s m j=0 w j g (αi) n (t j ) y ̟ (t m-j ), (22)
where g

(αi)
n (t j ) := [D αi t g n (t)] t=tj . The integrals given in ( 13) and ( 15) can be approximated in a similar way.

By writing y ̟ (t j ) = y(t j ) + ̟(t j ), we get:

Y ̟ N,m (n, i + 1) = Y N,m (n, i + 1) + e ̟ N,m (n, i + 1), (23)
where

Y N,m (n, i + 1) = T s m j=0 w j g (αi) n (t j ) y(t m-j ), (24)
e ̟ N,m (n, i + 1) = T s m j=0 w j g (αi) n (t j ) ̟(t m-j ). (25)
Thus the integral Y ̟ N (n, i + 1) is corrupted by two sources of errors:

• the numerical error which comes from a numerical integration method, • the noise error contribution e ̟ N,m (n, i + 1). Consequently, the estimation errors for the estimators given in Proposition 2 come from both the numerical errors and the noise error contributions in the discrete noisy case.

It is well known that if the value of T is set, then when T s tends to 0, i.e. m → +∞, the numerical errors tend to 0. In the next subsection, we are going to study the effect of the sampling period on the noise error contributions.

C. Influence of sampling period on noise error contributions

In this subsection, we consider a family of noises which are stochastic processes satisfying the following conditions:

(C 1 ) : for any t, s ∈ I, t = s, ̟(t) and ̟(s) are independent;

(C 2 ) : the mean value function of ̟(•) denoted by E[•]
belongs to L(I);

(C 3 ) : the variance function of ̟(•) denoted by Var[•] is bounded on I. Note that the white Gaussian noise and the Poisson noise satisfy these conditions. Then, we can give the following proposition.

Proposition 4: Let ̟(•) be a stochastic process satisfying conditions (C 1 ) -(C 3 ), and ̟(t j ), for j = 0, • • • , m, be a sequence of {̟(•)} with an equidistant sampling period T s . If D αi t g n (t) ∈ L 2 (I), then we have the following convergence in mean square of the noise error contribution in the integral Y ̟ N (n, i+1) given in [START_REF] Podlubny | Fractional Differential Equations[END_REF], for n = 1, . . . , N and i = 0, . . . , L: e ̟ N,m (n, i + 1)

L 2 (I) = == =⇒ Ts→0 T 0 D αi t g n (t) E [̟(T -t)] dt, (26)
where e ̟ N,m (n, i + 1) is given in [START_REF] Perruquetti | An algebraic approach for Humane posture estimation in the sagital plane using accelerometer noisy signal[END_REF]. Moreover, if ∀ t ∈ I, E [̟(t)] = 0, then we have:

e ̟ N,m (n, i + 1) L 2 (I) = == =⇒ Ts→0 0. ( 27
)
The proof of the previous proposition can be obtained in a similar way to the one given in [START_REF] Liu | Error analysis of Jacobi derivative estimators for noisy signals[END_REF]. Moreover, a similar result was studied using the non-standard analysis in [START_REF] Fliess | Analyse non standard du bruit[END_REF].

Consequently, according to the previous proposition, the noise error contributions can be increasing with respect to the sampling period.

Finally, let us mention that solving the linear system given in Proposition 2 in noisy case is related to the matrix perturbation theory. The accuracy and the stability of the proposed estimators not only depend on the noise error contributions, but also depend on the condition number of the associated matrix. This condition number depends both on the input and the output of the fractional order system and on the used modulating functions. In general, we should choose the modulating functions that can give a small condition number. This study is out of the scope of this paper.

V. SIMULATION RESULTS

In order to illustrate the accuracy and robustness with respect to corrupting noises of the proposed estimators, we present some numerical results in this section.

Let us consider a fractional order system defined by the following fractional differential equation:

∀ t ∈ [0, 8], a 0 D α0 t y(t) + a 1 D α1 t y(t) + a 2 D α2 t y(t) = u(t), (28)
where a 0 = 3, a 1 = 2, a 2 = 1, α 0 = 0, α 1 = 0.5 and α 2 = 1.5. We assume that the output is y(t) = sin(3t) + 1.

Hence, the initial conditions of y are not equal to 0. Moreover, the expression of the input can be obtained using (2) and the following formula (see [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] p. 83):

D αi t sin(3t) = 3t 1-αi Γ(2 -α i ) 1 F 2 1; 1 2 (2 -α i ), 1 2 (3 -α i ); - 1 4 3 2 t 2 , (29)
where p F q (c 1 , . . . , c p ; d 1 , . . . , d q ; •) is the generalized hypergeometric function given in [START_REF] Miller | An Introduction to the Fractional Calculus and Fractional Differential Equations[END_REF] p. 303.

In our identification procedure, we use the following modulating functions, the fractional derivative of which are simple to calculate:

g n (t) = (T -t) 6+n t 6+N +1-n , (30)
where n = 1, 2, • • • , N with N = 13. The value of T is taken to be equal to the time where we estimate the parameters. Let us recall that this kind of functions has been obtained when the algebraic parametric estimation technique was applied to the parameter estimation for signals described by differential equations [START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF].

In the two following examples, we estimate the parameters a 0 , a 1 and a 2 using the noisy observation of y where the noise is a high frequency sinusoidal noise and a gaussian noise, respectively. Moreover, we apply the trapezoidal rule to numerically approximate the integrals in our estimators.

Example 1. In this example, we assume that y ̟ (t j ) = y(t j )+0.5 sin(10 3 t j ) with T s = 0.01. We can see this discrete noisy signal in Figure 1. In our identification procedure, we take T = t i for t i ∈ [1.8, 8]. The obtained estimations and the associated relative estimation errors are shown in Figure 2 and Figure 3. Hence, we can see that the proposed estimators can cope with a high frequency sinusoidal noise.

Example 2. In this example, we assume that y ̟ (t j ) = y(t j ) + σ̟(t j ), where T s = 0.01, ̟(t j ) is simulated from a zero-mean white Gaussian iid sequence, and σ ∈ R * + is adjusted in such a way that the signal-to-noise ratio is equal to SN R = 22dB. This noisy observation is shown in Figure 4. The obtained estimations and the associated relative estimation errors are given in Figure 5 and Figure 6. We can see that the proposed estimators are robust against a gaussian noise.

VI. CONCLUSIONS

In this paper, the modulating functions method has been generalized to the on-line identification problem of fractional order systems. Using this method, the unknown parameters have been estimated by solving a linear system of algebraic equations involving the input and the noisy output. Thanks to the properties of modulating functions, we do not need to estimate the fractional derivatives of the output, to calculate the ones of the input. We do not need to know the initial conditions either. Moreover, the integral given in the estimators can reduce the noise effects due to a high frequency sinusoidal noise or a class of stochastic processes. The efficiency and the robustness against corrupting noises have been confirmed by numerical examples. In order to improve the robustness against noises, some methods such as the instrumental variable method will be applied [START_REF]Identification of continuous-time models from sampled data[END_REF]. It was mentioned that the proposed estimators also depend on the choice of modulating functions. This problem will be studied in the future work. Moreover, the estimation of the fractional derivative orders in a time-delayed fractional order system will be considered.

Fig. 1 .

 1 Fig. 1. The noise-free output and its noisy observation.

Fig. 2 .

 2 Fig.2. The exact parameters and their estimations with a sinusoidal noise.

Fig. 3 .

 3 Fig. 3. Relative estimation errors in the sinusoidal noise case.

Fig. 4 .

 4 Fig. 4. The noise-free output and its noisy observation.

2 Fig. 5 .

 25 Fig.5. The exact parameters and their estimations with a gaussian noise.

Fig. 6 .

 6 Fig. 6. Relative estimation errors in the gaussian noise case.

More generally, the noise is a stochastic process, which is integrable in the sense of convergence in mean square[START_REF] Liu | Parameters estimation of a noisy sinusoidal signal with time-varying amplitude[END_REF].


Ahmed Gamal Eldin

Guillaume Charpiat

Xavier Descombes

Josiane Zerubia

Ahmed Gamal-Eldin

An efficient optimizer for

Keywords: Point process, multiple birth and cut, non-convex optimization, multiple object detection

come

INTRODUCTION

In the last decade, multiple object detection has evolved in both models and optimizers. We find two different mainstreams from the aspect of object modeling. Inverse methods follow a bottom-up approach: they construct the objects from primitives, like corners, edges, feature points or segmented object parts, by a rule based model or a previously trained classifier. [START_REF] Pham | Fast polygonal integration and its application in extending haar-like features to improve object detection[END_REF] Recognition with these solutions can be notably quick due to their sequential workflow, however, their performance may drastically decrease if the features cannot be reliably detected. On the other hand, direct methods assigning a fitness value to each possible configuration, and an optimization process attempts to find the configuration with the highest confidence. This way, flexible object appearance models can be adopted, and it is also straightforward to incorporate prior shape information and object interactions into the model. However, the optimization process can be a bottleneck for many sophisticated models. The richer the model is, the more complex its optimization becomes. This phenomenon makes most of the research focus on very simple models that can easily be optimized, even if these models can be over-simplistic for some applications. Markov random field (MRF) models are life examples that suffer this type of problem. While higher interactions with sophisticated potential terms exhibit more general modeling capabilities, their optimization is very hard. Consequently, most of the popularly used models are simple models that can be efficiently and rapidly optimized. But these models can be oversimplistic for some applications. Considering MRFs, the most widely used class is of first order with limited set of potential terms, because this class can be efficiently optimized with graph cut algorithms. [START_REF] Boykov | Markov random fields with efficient approximations[END_REF][START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF] More sophisticated models require stochastic samplers such as Gibbs and Monte Carlo Markov Chain (MCMC) algorithms. In this work we consider a multiple object detection framework known as Marked point process. This framework has shown a great potential in remote sensing applications, such as building detection, tree detection, road detection, and many others. The strength of this model lies in its capacity to incorporate a wide variety of data terms and also a wide range of prior information. Data terms for object detection can range from simple correlation to using an SVM classifier trained on SIFT features for specific object category. The prior term also can range from simple non-overlapping condition to prior information about minimal distances and angles between objects.

Here the optimization issue plays a crucial role due to the high dimension of the configuration space. Until very recently, the only existing methods to optimize MPP models were fully stochastic. The initial optimizer is Birth-and-Death, [START_REF] Geyer | Simulation Procedures and Likelihood Inference for Spatial Point Processes[END_REF] then came the Reversible Jump Markov Chain Monte Carlo (RJMCMC) [START_REF] Green | Reversible Jump Markov Chain Monte Carlo computation and Bayesian model determination[END_REF] and recently came the Multiple-Birth-and-Death (MBD). [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF] While each of these optimizers represents on its own a major advance in both the optimization world and in the application world, these algorithms are fully stochastic and are embedded in the simulated annealing scheme. Although those algorithms may be able to optimize very complex models with huge configuration space, they are very slow in practice.

A recent work [START_REF] Gamal Eldin | Multiple birth and cut algorithm for point process optimization[END_REF] was the first to show the possibility of optimizing such models (MPP) with a mixed semideterministic optimizer using graph-cut. This new algorithm is known as Multiple Birth and Cut (MBC). This method holds two parts, the stochastic part to explore a very large configuration space, and the deterministic part which makes the optimal selection between an existing configuration and a newly proposed one.

In this paper, after a presentation of our Point Process Model, we propose methods to boost the efficiency of the MBC optimizer. First, we propose a simpler and faster way to speed up the birth step of this algorithm than in. [START_REF] Gamal Eldin | A fast multiple birth and cut algorithm using belief propagation[END_REF] Next, we introduce concepts from the stochastic MCMC optimizer, in a simpler way that will again improve the speed of the optimizer. Then in section 5, we study theoretically the complexity to express convergence properties, and present how this algorithm scales on a synthetic problem, compared to the existing algorithms. Finally, in section 6, we show some of the many possible applications to which this type of modeling can be applied.

POINT PROCESS MODEL

We consider that point process models in image processing can be divided into two classes. The first class covers models with simple interactions. Examples for this class could deal with trees, flamingos, cells, and many others. The second class is appropriate for modeling scenes of objects with complex interactions (priors), such as roads. For example, road segments have to be connected, and the intersection angle can not be less than 5 • . In this paper we focus on the first class. In the following, we briefly introduce point process models.

We consider a point process X living on a space K. K is the image and is defined by

K = [0, I max ]×[0, J max ],
it is a closed connected subset of R 2 . The point process X is a mapping from a probability space (Ω, A, P) to the set of configurations of points on K. The points on their own represent the objects' centers, while the marks hold the geometrical parameters of these objects. Each object ω i consists of a point x i and a mark m i associated to this point; let M is the mark space. The configuration space is then defined as:

Ω = ∞ n=0 Ω n , Ω n = {{ω 1 , . . . , ω n } ⊂ K × M } , (1)
where Ω n is the subset of configurations containing exactly n objects. The number n of objects to be detected is a random variable which is unknown.

The density of the process is defined by a Gibbs density on the configuration space. The Gibbs model provides a powerful tool for modeling a scene of objects and interactions between them. The energy is written as the sum of potentials over interacting objects (cliques). The density is therefore defined by:

f (ω) = 1 Z exp[-U (ω)] = 1 Z exp[-(U d (ω) + γ p U p (ω))] , (2)
where U p is the prior energy which takes into account the interactions between geometric objects, U d is the data energy to fit the configuration to the image, Z is the partition function and γ p is the weight of the prior term.

Data Term

As we mentioned in the introduction, the MPP framework does not impose any limitations on the selected method for evaluating the fitness of a proposed object. Given the independence of the data term of each object, the data term energy of a configuration ω is given by U d (ω) = ωi∈ω u d (ω i ) The term u d (ω i ) evaluates from a data point of view the relevance of the proposed object ω i . The object contains information on both its location and its shape.

Prior Term

The capacity to incorporate prior information (for regularization) is one of the major advantages of MPP models.

Since in this work we only consider the first class of MPP models, we have a simple prior term which consists in penalizing overlapping objects, as objects in the selected applications in the real world do not overlap. Let A(ω i , ω j ) ∈ [0, 1] be the overlapping coefficient between two objects. As required by the MBC algorithm for the optimization of the first class of MPP models, the interaction cost should be zero or infinity. The prior energy of a local configuration {ω i , ω j }, where ω i ∼ ω j ( ∼ indicates neighbor objects) is given by:

u p (ω i , ω j ) = 0 if A(ω i , ω j ) < 0.1 ∞ if A(ω i , ω j ) ≥ 0.1, (3)
which means that we tolerate small overlapping up to 10%, otherwise it costs infinity. The total prior configuration energy is given by U p (ω) =

(ωi∼ωj )∈ω u p (ω i , ω j ).

OPTIMIZATION

Now that we have defined the density by an energy function U (.), we want to compute the Maximum Likelihood (ML) estimate with respect to the density f . The ML estimate corresponds to the configuration that matches best our model and the input image. The estimate of the global mode of a density is given by:

ω = argmax ω∈Ω f (ω) = argmin ω∈Ω U (ω)
Finding the minimum (or maximum) for this highly non-convex function (equation 2) with many local minima requires a global optimization method.

What makes a good optimizer?

In this section we will discuss what we consider to be the main characteristics of optimizers 9, 10 for non-convex functions in the context of MPP models.

1. Global versus local proposals: Local proposal * means exploring the fine details of the modes, while global proposal means exploring the different modes. For a multimodal distribution, and with only local proposals, based on the initialization (starting point), only one mode will be discovered, while using only global proposals will not capture the details the modes (it will be very slow to capture the detailes). For either simulation or optimization, given that the considered density is multimodal, both these properties are required.

Simple versus multiple perturbations:

The perturbations are simple or multiple based on how the configuration changes between two iterations. If only one object is modified (changed, removed or added), then it is a simple perturbation. If many objects are modified per iteration, then it is considered as a multiple perturbation.

Holding those characteristics is essential to make a fast optimizer. It requires global and local perturbation kernels, and also to be able to make multiple perturbations (simple perturbation is implicite).

Kernels Evolution in Point Process models

• In a multiple object detection problem, the dimensionality of a set is unknown in advance, as it is the number of objects to be detected. This means that only samplers that can make dimensional jump can be used. This was only possible after the development of a birth-death kernel by Geyer et al. [START_REF] Geyer | Simulation Procedures and Likelihood Inference for Spatial Point Processes[END_REF] • MCMC is a flexible sampler, it can incorporate different kernels that are problem specific, but it can not make a dimensional jump. By integrating the birth-death kernel into MCMC, Green obtained a new sampler that can also make dimensional jump. This sampler is known as RJMCMC.

• Using the previous samplers, we can only make simple perturbations. The problem is that in practice those samplers are very slow. Due to the work of Descombes et al., [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF] it became possible to make multiple perturbations. They developed a multiple-birth-death kernel, and this is global multiple perturbations kernel.

• All the above samplers are fully stochastic embedded in a simulated annealing scheme. It is only very recently that a semi-deterministic optimizer was introduced by Gamal-Eldin et al. [START_REF] Gamal Eldin | Multiple birth and cut algorithm for point process optimization[END_REF] This optimizer is known as Multiple Birth and Cut (MBC). It also has a multiple-birth-death kernel for global multiple perturbations.

NEW MBC ALGORITHM

We start from the MBC algorithm presented in. [START_REF] Gamal Eldin | Multiple birth and cut algorithm for point process optimization[END_REF] We propose a more efficient multiple-birth-death kernel than the one proposed in. [START_REF] Gamal Eldin | A fast multiple birth and cut algorithm using belief propagation[END_REF] The proposed kernel makes more efficient global perturbations, and explore more efficiently the configuration space, while being very simple. We next show how local perturbation kernels can be integrated in the MBC algorithm to make it more efficient.

New selection method

We start by analyzing the behavior of the selection method based on belief propagation presented in. [START_REF] Gamal Eldin | A fast multiple birth and cut algorithm using belief propagation[END_REF] We consider a 1D problem as illustrated in figure 1(a). Let ω [n] be the current configuration, presented in blue in figure 1. Assume that the two empty spots should contain objects. Over each of the two empty spots a packet of objects is proposed. The packets are presented in red in the figure. The aim is to select the best candidate object inside each packet X i . There exist some interactions (presented by the gray dotted lines) between some objects of different packets. Each object has a data term (cost) and a prior term (interaction cost). Let the set of proposed packets be defined by Γ = {X 0 , X 1 , . . . , X n }. The belief propagation algorithm makes the best selection from Γ (global optimal), based on the data terms and the prior terms. The algorithm ends by selecting the two candidate objects {ω 1 1 , ω 1 2 }, shown in green in figure 1(b). While the MAP estimate using belief propagation gave the optimal selection from Γ, those selected objects were not necessarily optimal from the point of view of the current configuration ω [n] . As figure 1(c) illustrates, what was optimal for the current configuration ω

[n] was {ω 2 1 , ω 3 2 } and not {ω 1 1 , ω 1 2 }. The difference is, {ω 2 1 , ω 3 2 } is locally optimal, while {ω 1 1 , ω 1 2 } is optimal over Γ.
The aim is to minimize the total energy. If we consider only the data terms,

u d (ω 1 1 )+u d (ω 1 2 ) < u d (ω 2 1 )+u d (ω 3
2 ), so ω 1 1 and ω 1 2 are perfect candidates, while ω 2 1 and ω 3 2 are good candidates. If we consider the total energy, prior term and data term,

u d (ω 1 1 ) + u d (ω 1 2 ) + u p (ω 1 1 , ω 1 2 ) > u d (ω 2 1 ) + u d (ω 3 2 ).
We propose to make the selection only locally, and remove objects that violate the non-overlapping condition. We re-propose new objects instead of the rejected ones. It costs much more to get a perfect candidate than to get a good one. So keeping one of the two perfect candidates and re-proposing instead of the other is more efficient. In real applications, the percentage of rejected objects ranges between 5% and 10%. We keep iterating the algorithm until re-proposing this small percentage. The main reasons for being faster are:

1. It takes more time to propose an excellent candidate than a good one.

2. Let the size of Γ be m × n, where m is the number of packets and n is the number of objects per packet.

The complexity of belief propagation is O(mn 2 ), while for local selection it is O(mn).

3. Proposed configurations can be very dense, and they can form loops.

4. The belief propagation algorithm require extra computational cost to ensure that the proposed packets do not form any loop.

Local Perturbations Kernels

The existing multiple perturbation optimizers hold a single kernel, a birth kernel that proposes new objects.

What is missed here is making local proposals.

Consider object ω i to be the best candidate we obtained from the birth kernel to detect object o i . Object o i has a position x i and a set of marks m i . If we consider the ellipse model in R 2 , the set of parameters of ω i is Θ = {x i , y i , a i , b i , ρ i }. Given that object ω i is a good candidate for the detection of object o i , this means that Θ oi is close to Θ ωi . It is much easier to go from Θ oi to Θ ωi (or closer) than to propose a new object that becomes a better candidate for object o i . This can be obtained via perturbation kernels.

Let us consider four perturbation kernels. One for perturbing the object position, a second to perturb the major axis, a third for the minor axis and a fourth one for the angle. From the current configuration ω [n] we randomly select a subset ω of objects. On each object ω i ∈ ω we apply a perturbation kernel, ω perturbations ---------→ ω . More than one kernel can be applied to each object. Kernels are selected randomly. Giving different weights to the kernels may be beneficial. This type of perturbation is usually referred to as local random walk. Those kernels were originally proposed in the context of RJMCMC sampler. [START_REF] Perrin | Forest resource assessment using stochastic geometry[END_REF] While algorithm is not and MCMC algorithm, but we were able to take advantage of the many kernels idea. The new MBC algorithm is summarized in algorithm 1.

Algorithm 1 New Multiple Birth and Cut

1: n ← 0 , R ← const1 2: generate ω , ω [0] ← ω 3: repeat 4:
Birth:

5:

Sample u ∼ U (0,1)

6:

Based on u, select a kernel randomly (including the birth kernel)

7: ω ← Apply selected kernel 8: ω ← ω [n] ∪ ω 9: Cut: ω [n+1] ← Cut(ω [n]
∪ ω ) (optimize with graph cuts) 10: until converged Global proposals are mostly needed at the starting of the algorithm. We force for the first few iterations to only use of the birth kernel. Thereafter we alternate between all the kernels. The proper selection of the models parameters and dynamic kernels selection are based on the following analysis.

ALGORITHM ANALYSIS

5.1 Theoretical analysis 5.1.1 Notations.

Algorithm 1 alternates randomly two kinds of steps: global exploration, with probability p G , and local exploration, with probability p L = 1 -p G . A global exploration consists in picking randomly α G N locations in the image (using the birth map, where N is a rough estimate of the number of objects to be found, and α G ∈ [0, 1] a proportion) and in testing n G random marks at each of these locations. A local exploration consists in picking randomly α L N already detected objects, and in testing n L random small variations of their marks (i.e. exploring direct neighbors in the discrete space of marks M of angle, axes' lengths and precise location).

Preliminary step.

Before applying Algorithm 1, a series of successive global explorations is performed, in order to ensure that most of the image domain K has been covered and that the number of objects detected (even if with wrong marks) is close to N . The average proportion of the image still uncovered after k such steps, i.e. the average proportion of regions with area |K| N not visited yet, is (1 -α G ) k . For α G = 60% e.g., k = 14 steps are sufficient to ensure that 99.9% of the image has been explored.

Global vs. local explorations.

Now, at step s of Algorithm 1, the average number of random tests already performed in the neighborhood (of

area |K| N ) of any image point x is v G s with v G := α G p G n G ,

and the average number of local explorations around any detected object is similarly

v L s with v L = α L p L n L .

Thus the ratio of local tests over global tests at any location is v

L v G = 1-p G p G α L n L α G n G .
However, the probability that a global test is successful decreases with time: by construction it is at most the probability that a random try in M performs better than the previous v G s -1 ones, which is 1 v G s . Consequently, the expected number of local explorations between two big moves in M increases with time and is v L v G 2 s , which means that local minima are more thoroughly explored with time.

Time complexity and Percolation.

We will now study the time spent at each step, in order to optimize the algorithm w.r.t. parameters v and α. The cost of any exploration step, global or local, is of the form αnN , plus the cost of applying graph cut to a graph of approx. 2αN nodes. The expected complexity of graph cut is much lower than its worst-case complexity, and is usually regarded as between linear and quadratic in the size of the graph, depending on problems. The graph here however is generally not connected and thus the complexity depends on the size of its connected components (the smaller, the lower). Since the graph is made of randomly selected objects with density α, we are precisely interested into percolation properties, to obtain the typical number and size of connected components. [START_REF] Kesten | Percolation Theory for Mathematicians[END_REF] The critical site-percolation density α c , above which arbitrarily big components arise, depends on the neighborhood size and on the type of lattice the graph is extracted from. In our application case, the number of neighbors is at most 6, and one can lay an hexagonal lattice on the image, with step size the typical object size. We could also consider the Voronoï diagram associated to the objects already detected, which is relatively similar. In this hexagonal setting, α c is known to be approx. 0.697... . Since α has to satisfy α < α c , a very cautious choice is to set α 60%. For such an α, the phase is subcritical and the probability of connected regions is known to decrease exponentially with their area A, as e -A/R , where R depends on α (smoothly for α far from α c ). The typical region size is then R, and, by computing A A 2 e -A/R dA / A e -A/R dA , we obtain that the average graph cut cost (supposing quadratic complexity) is proportional to R 2 . One can bound R by R M = R(0.6) if we restrict α 60%. Then the average total graph cut cost at step s, which involves 2αN objects in about αN R connected components, is bounded by αN R M . Consequently, the complexity of step s is bounded by C S = αN (n + R M ). Note that this bound does not depend on s. The expected clock time at the end of step s is thus

t = (p G C S (G) + p L C S (L)) s = (v G + v L + R M [p G α G + (1 -p G )α L ]) N s
which has to be compared to v G s and v L s, the average number of explorations in the neighborhood of any image point.

Estimated time and precision when in an attraction basin.

Let us consider the case of a locally simple problem, in the sense that objects are located away enough from each other, in order to prevent any mutual interaction. Under this hypothesis, any suitable data energy will satisfy that each individual potential admits a unique minimum, i.e. that the attraction basin of any ground-truth object includes all possible marks located in the neighborhood of this object. We can then estimate the time needed to find the global minimum in the mark space M with the method above. After x random global explorations in an image region of area |K| N , the minimum is at distance

t = (1 + R M n G )(v G s) + (1 + R M n L )(v L s) N , which is solved by v G = v L ,
and by high values of n G = n L >> R M to reduce the relative cost of graph cuts (we choose n = 10 in practice since R M observed is already low). The corresponding average time to find the optimal mark for an object is 2 |M |N . More precisely, the probability that the optimal mark is not found for this object

after time t is p N = |M |! (M -t/N )! M t/N ≈ e -
t = γ |M |N is (1 -e -γ 2 ) N ≈ e -N exp(-γ 2 )
; conversely the time needed to reach a given probability p is t = γ(p)N |M | with γ(p) = -ln -ln p N . Furthermore, the expectancy of the proportion of the objects not found yet after time t = γ |M |N is e -γ 2 (in the case of a spatially uniform birth map). Thus a reasonable total time to spend in the locally convex case to expect all but one objects and their optimal marks is then |M |N √ ln N .

Fixing local minima.

The results above remain valid for any problem provided objects are correctly localized from the first steps. However, it may happen, in the case of an object category with very-varied sizes, that two contiguous small objects are proposed instead of a big one. In such a case, one needs a global exploration at the correct location with reasonably good marks. The required quality of the birth depends on the relative energy cost of the wrong couple of small objects w.r.t. the one of the correct detection. Given a lower bound B > 0 on this difference (based on misplaced edges e.g.), it is sufficient that a global exploration finds a mark with energy less than B + the optimal one. Let p GM be the proportion of such marks within the local mark space M . The average time to find one of them is 1 2p GM N . With similar techniques as before, the time required to fix, with probability p, simultaneously k such local minima occurring in the image, is ln(1 -p 1/k ) 1 2p GM N . Note that another, faster possibility would be to add a merge-and-split kernel dedicated to this kind of minima.

Semi-local minima.

Depending on the energy minimized to detect objects, and in the case of overlapping neighborhoods of objects, bad initial explorations may result in wrong detections half-way between two real objects. Such wrong detections are fixed quickly by Algorithm 1, when a global exploration is performed at one of these locations, which happens as frequently as 2p G α . What is more problematic is when such wrong detections form chains of misplaced objects, each one preventing its neighbors from evolving to its right place. If the chain is a tree of diameter d, then the expected time to remove it is only d p G α . However if the chain comprises a cycle of length c, the only way to break it (if this local minimum is strong) is to wait for a simultaneous global exploration of all c sites, which happens as rarely as p G α c . Nevertheless, the appearance of such a cycle is possible only during the very few first steps, and the probability of such an event is as low as p c α c where p is the probability of one misplacement. Moreover, whenever a chain is broken, it cannot form again. Due to space limitation, some of the simple mathematical details were skipped.

How the algorithm scales

A comparison now is essential to see how the proposed kernels affect the speed of convergence. The comparison will concern the different versions of the MBC algorithm. For the seek of completeness, we also add the MBD algorithm which is the fastest one in the stochastic family for this class of MPP models. [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF] For details of the MBD algorithm, please refer to. [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF][START_REF] Descamps | Automatic flamingo detection using a multiple and death process[END_REF] Other models for object counting exist, but we did not compare with those methods for two reasons: (1) the focus of this paper in on the optimization of MPP models, (2) methods such as 14 requires a training phase while our does not require.

We tested the four algorithms on three samples containing respectively 300, 1000 and 10000 objects. These samples are generated by the same ellipse model introduced in. [START_REF] Gamal Eldin | Multiple birth and cut algorithm for point process optimization[END_REF] The 300 object sample is shown in figure 2(a), and the detection using the new MBC algorithm is shown in figure 2(b). The comparison will present the evolution of the total energy of configurations versus time for each algorithm, as well as the object detection rate. In figure 3 we present the result on the largest sample of 10000 objects. In this figure we refere to the naive MBC by mbc1, mbc2 referes to MBC using belief propagation, mbc3a is our algorithm using only the new birth kernel and mbc3b refere to our algorithm using all kernels.

We use the following set of parameters: for the MBD algorithm, temperature T = 1/50, ∆T = 0.9997, ∆β = 0.999. For the MBC algorithm, the number of packets (objects for the naive MBC) is R = 2000, while for the first few iterations, we used R = 3 × R (for faster convergence), the packet size is 8 (no packets in the naive MBC).

For this configuration, figure 3(a) shows how the energy of configurations evolving with time for the four algorithms. The true value of the energy is E(ω) = -7740.95. From the energy curves shown in figure 3(a,b), we can conclude the following:

• MBD converges faster that the naive MBC, but setting its parameters is not straightforward.

• Naive MBC (mbc1) will converge to the correct minimum, but it is slow.

• MBC using belief propagation (mbc2), is faster than the MBD and the naive MBC, but after a while its converges slowly.

• By adding only the new birth kernel, the algorithm (mbc3a) is much faster than the previous ones.

• By including all the proposed kernels, the algorithm (mbc3b) becomes very fast. Compared to the MBC using belief propagation, our algorithm is around 20 time faster.

From the detection rate curve shown in figure 3(b), we conclude that: • All the algorithms find the correct number of objects in around 10000 sec.

• Our new algorithm using all kernels has the least variance to the true number of objects. In 300 sec, with a positive error of 2%, and then starts decreasing again to the correct value.

• The detection speed of the MBD algorithm is quite good after 1000 sec.

We can get a very good result from our new algorithm just after 100 sec for 10000 objects, while for the MBD we have to wait till 1000 sec.

Algorithm Parallelization

Since in this paper we consider only the first class of MPP models, happily it can be very easily parallelized. Given that the objects we are modeling, for e.g. flamingos, are independent, the resultant configuration from the optimization algorithm also should contain independent (non-overlapping) objects.

Given an input image, it should be split into n partitions, where n can be the number of cores in a processor, or number of machines used to solve the problem. Two issues to consider for splitting this image:

1. Each time we divide the image, there should be an overlap equal to the size of the largest existing object.

2. In many real cases, dividing the image into two equal halves at each split may not be the optimal choice for load balance. We propose to first calculate the birthmap, then make the split based on this map.

After splitting, each part is processed independently. The last step of the algorithm will be to join the result on the different parts. We propose to make the optimal join using graph cut, as it is done in every iteration in the Cuts step of the MBC algorithm.

To illustrate, consider the small sample of a flamingo colony shown in figure 4(a). This image is split vertically into two partitions, and each is optimized independently using the MBC algorithm. In figure 4(b), we overlap the detection result of the two partitions, one in pink and the other one in green. The common margin contains two detection result. In this example, nine flamingos are detected by both runs. We construct a graph, and to apply the same idea used in any of the MBC algorithms to select between two sets of configurations. This way of parallelization has the following advantages:

1. The cost of the fusion step is the of cost optimizing one small graph. 4. If an optimization algorithm other than MBC is needed, it can be used without any change, the fusion part will be the same.

APPLICATIONS

Recent technological advances in the development of airborne sensors have significantly increased the number of available satellites as well as their spatial and spectral resolutions. These advances, combined with the need for an effective environmental control, have favored the emergence of new remote sensing applications. Here we present results on three different applications.

Ecological example 1: Flamingo counting

The unprecedented biotic crisis, due to human activities, follows from the destruction and fragmentation of natural habitats, overfishing, pollution and climate change. We consequently need long term monitoring systems allowing to estimate each year the number of breeding flamingos in the colonies and their breeding success. The development of an automatic counting system by photographic analysis of flamingo colonies could be a powerful and economic tool. We apply our MBC algorithm using the ellipse model developed in [START_REF] Descamps | Automatic flamingo detection using a multiple and death process[END_REF] to accomplish this automated counting. An aerial image of a small flamingo colony taken in Turkey is shown in figure 5.

Ecological example 2 : Tree counting

Brazilian eucalyptus plantations are among the most productive ecosystems in the world, providing wood sources for industry and energy. Monitoring these plantations on large surfaces is rather difficult and requires a lot of manpower. The aim is detecting individual tree crowns in localization and in size on large areas automatically, and provides important statistics for planting work.

In figure 6 we present a high resolution satellite image of eucalyptus plantation. We also considered the ellipse model while setting the acceptable overlapping coefficient to zero to obtain the result shown in figure 6(b).

Surveillance: Boat detection

Here we present a more complicated example with a peaky energy. Given a satellite image of a port presented by figure 7(a), the aim is to detect and extract the boats. This energy is peaky because of objects' orientation: there exists a special structure (alignment) in the data. We use ellipses to approximate boats' shapes. In figure 7(b), we present the detection result, assuming uniform prior on the orientation of the ellipses θ ∼ U [0,π[ . The detection result is not that good, because of the complexity of finding this global minimum. One possible solution, is to run the algorithm for much longer time. Another solution will be to reduce the complexity of this problem by using a non-uniform prior on the orientation parameter. Another alternative would be to add an extra term to the prior term for alignment between objects. Using a non-uniform prior on the orientation parameter, we get the result shown in figure 7(c). Figure 7(d) shows the detection result by adding an alignment term to the prior energy. [START_REF] Ben Hadj | Parameter estimation for a marked point process within a framework of multidimensional shape extraction from remote sensing images[END_REF]

CONCLUSION AND FUTURE WORK

In this paper, we discussed the essential characteristics that should be considered to design efficient optimizers, while focusing on a selected class of models, namely Point Process Models. We presented the limitations of existing birth-and-death kernels, and, to overcome them, we introduced new kernels, both local and global.

In a theoretical analysis, we demonstrated the importance of both global and local perturbations. This analysis also showed how to select optimally the optimizer's parameters; in particular it quantified the optimal times spent in local and global explorations, in order to obtain the fastest convergence rates. We validated our algorithm on three different real applications, and showed its efficiency, as well as how it scales with the problem's size on synthetic data.

The MBC optimizer has proven to be simple and modular, while being very efficient. The proposed kernels are very simple and can be implemented easily. Moreover, adding new kernels to adapt the algorithm to specific applications, or even to tune the current algorithm, is possible without a lot of burden. While having all those simply designed kernels, the Cut step which is based on a binary graph cut, guarantees a global optimal selection. 2, 7 Future work will consider extending this algorithm to solve more complex models, such as Point Process Models with sophisticated interactions.

Figure 1 :

 1 Figure 1: (a) A synthetic 1D problem. Current configuration ω [n] shown in blue. Candidate objects shown in red. Interaction terms are presented with dotted gray lines. (b) Best candidates globally are shown in green. (c) Best candidates locally are shown in green.

Figure 2 :

 2 Figure 2: (a) Sample of 300 objects using the ellipse model after adding a Gaussian noise. (b) The detection result using the new MBC algorithm.

Figure 3 :

 3 Figure 3: (a) Shows the energy evolution versus time on the 10000 objects' configuration. (b) Shows the number of objects detected as a function of the time spent (in log scale).

Figure 4 :

 4 Figure 4: (a) A small part of an aerial image of a flamingo colony. (b) Synthetic result on this sample, each color represent objects obtained independently on two cores.

Figure 5 :

 5 Figure 5: (a) An aerial image of a small flamingo colony taken in Turkey in 2006 (b) The detection result, each flamingo is surrounded by a pink ellipse.

Figure 6 :

 6 Figure 6: (a) Brazilian eucalyptus plantation (b) The detection result, each tree is surrounded by a pink ellipse.

Figure 7 :

 7 Figure 7: (a) photograph of vessels in France c CNES. (b) Detection result when using a uniform prior on boats orientations. (c) Detection result when using a non-uniform prior on the orientations. (d) Result when adding an alignment term to the prior energy.

  |M | x on average and thus then reachable in |M | x local explorations. Minimizing x + |M | x leads to x = |M |. Thus we aim at minimizing t such that v G s |M | and v L s |M |, w.r.t. all parameters. This leads to minimizing


J Martinod

D Hatzfeld

J P Brun

P Davy

P Gautier

Continental collision, gravity spreading, and kinematics of Aegea and Anatolia

We have carried out experiments using a layered medium of sand and silicone to investigate the lateral extrusion of a material which spreads over its own weight while being compressed by the advance of a rigid indenter. Boundary conditions in the box mimic those prevailing in the Anatolian-Aegean system. Both shortening in front of the rigid piston, which models the northward motion of Arabia, and extension resulting from the gravity spreading of the sand-silicone layer are necessary to initiate the lateral extrusion. Strike-slip faults accommodate the lateral escape and link the normal faults accompanying gravity spreading with the thrust faults in front of the rigid indenter. Strike-slip faults begin to accommodate extrusion at a late stage in the experiments after the normal and thrust faults have developed. Experiments also show that the initial geometry of the boundary of the spreading layer may result in the formation of two arcs behind which material extends, in a manner analogous to the Hellenic and Cypriot arcs, without invoking a rheological change at the junction of the two arcs. The experiments also suggest that southward motion of the eastern part of the spreading region is compensated by the northward advance of the piston, which is a possible explanation for the slower movement of the Cypriot arc compared to the Aegean arc.

Introduction

While the present-day kinematics of the Anatolian-Aegean system are now well constrained by geodetic Global Positionning System (GPS) and Satellite Laser Ranging (SLR) measurements [e.g., [START_REF] Noomen | Crustal deformation in the Mediterranean area computed from space-geodetic observations[END_REF][START_REF] Oral | Global positioning system (GPS) measurements in Turkey (1988-1992) : Kinematics of the Africa-Arabia-Eurasia plate collision zone[END_REF][START_REF] Veis | Tectonic strain in Greece from geodetic measurements[END_REF][START_REF] Le Pichon | Geodetic determination of the kinematics of central Greece with respect to Europe: Implications for eastern Mediterranean tectonics[END_REF][START_REF] Reilinger | Global positioning system measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone[END_REF] (Figure 1), the tectonic processes responsible for the westward motion of the Anatolian-Aegean plate are still a subject of debate. Following [START_REF] Mckenzie | Active tectonics of the Mediterranean region[END_REF], most authors consider the westward movement of Anatolia with respect to Eurasia to be the result of the continental collision between Eurasia and Arabia, because extrusion minimizes the thickening of the continental crust in front of the Arabian indenter. It is worth noting, however, that if the Anatolian-Aegean plate bears a compressive tectonic regime on its eastern side resulting from the Arabia-Eurasia continental collision, it also bears extension deformation in its western half. Although the present-day rapid southwestward motion of the Aegean with respect to stable Europe (Figure 1) is consistent with the lateral escape of Anatolia to the west [e.g., Le Pichon et al., 1995;[START_REF] Armijo | Quaternary evolution of the Corinth Rift and its implications of the late Cenozoic evolution of the Aegean[END_REF], the thinning of the Aegean plate developed in response to local forces, such as the forces generated by the subducted Mediterranean oceanic lithosphere or those resulting from the topography and thickness of the Aegean plate at the end of the Hellenic orogenesis [e.g., [START_REF] Berckhemer | Some aspects of the evolution of marginal seas deduced from observations in the Aegean region[END_REF][START_REF] Le Pichon | The Hellenic arc and trench system: A key to the neotectonic evolution of the eastern Mediterranean area[END_REF][START_REF] England | Length scales for continental deformation in convergent, divergent, and strike-slip environments: Analytical and approximate solutions for a thin viscous sheet model[END_REF][START_REF] Sonder | Effects of a temperature-dependent rheology on large-scale continental extension[END_REF][START_REF] Meijer | Temporal variation in the stress field of the Aegean region[END_REF][START_REF] Lundgren | A geodynamic framework for eastern Mediterranean kinematics[END_REF]]. This extension first developed independently of the Arabia-Eurasia collision. As a matter of fact, it started in Cyclades as early as early Miocene [START_REF] Lister | Metamorphic core complexes of Cordilleran type in the Cyclades, Aegean Sea, Greece[END_REF][START_REF] Jolivet | 3-D kinematics of extension in the Aegean from early Miocene to present: Insight from ductile crust[END_REF][START_REF] Gautier | Ductile crust exhumation and extensional detachments in the central Aegean (Cyclades and Evvia islands)[END_REF], which is more than 10 Myr older than the Arabia-Eurasia collision [START_REF] Gautier | Timing, kinematics, and cause of Aegean extension: A scenario based on a comparison with simple analogue experiments[END_REF]. The arcuate form of the Hellenic trench turning from NNW-SSE in western Greece to NE-SW off Rhodos, the Neogene clockwise rotations in western Greece, the counterclockwise rotations in western Turkey [START_REF] Kissel | The Tertiary geodynamical evolution of the Aegean arc: A paleomagnetic reconstruction[END_REF], and the present-day extension parallel to the trench in the Hellenic arc [e.g., [START_REF] Hatzfeld | The strain pattern in the western Hellenic arc deduced from a microearthquake survey[END_REF][START_REF] Hatzfeld | Microearthquake seismicity and fault plane solutions in the southern Aegean and its tectonic implications[END_REF][START_REF] Kahle | The strain rate field in the eastern Mediterranean region, estimated by repeated GPS measurements[END_REF]] highlight the influence of gravity spreading of the continental plate above the eastern Mediterranean in the internal deformation of the Aegean [e.g., Le Pichon and [START_REF] Le Pichon | The Hellenic arc and trench system: A key to the neotectonic evolution of the eastern Mediterranean area[END_REF][START_REF] Sonder | Effects of a temperature-dependent rheology on large-scale continental extension[END_REF][START_REF] Hatzfeld | Gravity spreading: A possible contribution to the geodynamics of the Aegean[END_REF][START_REF] Gautier | Timing, kinematics, and cause of Aegean extension: A scenario based on a comparison with simple analogue experiments[END_REF].

Figure 1. Tectonics of the Anatolian-Aegean area. Present-day displacements are obtained

from [START_REF] Reilinger | Global positioning system measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone[END_REF]. Rotations are from the paleomagnetic data presented by [START_REF] Abrahamsen | Palaeomagnetic timing of the rotation and translation of Cyprus[END_REF], [START_REF] Kissel | The Tertiary geodynamical evolution of the Aegean arc: A paleomagnetic reconstruction[END_REF], [START_REF] Kissel | First paleomagnetic evidence for a post-Eocene clockwise rotation of the western Taurides thrust belt east of the Isparta reentrant (southwestern Turkey)[END_REF], and [START_REF] Speranza | Paleomagnetic evidence for Cenozoic clockwise rotation of the external Albanides[END_REF].

In the following, we present simple analogue experiments that simulate the Arabia-Eurasia collision, the Aegean extension, and the lateral escape of Anatolia, in order to study two points:

1. The first point is the effect of the compression applied on the eastern border of Anatolia, and of the forces generated near the Aegean that favor extension, on the motion of Anatolia. Experiments suggest that the lateral escape of Anatolia results from both the Arabia-Eurasia collision and the Aegean back-arc tension forces and that lateral escape could not occur without the two phenomena acting together.

2. The second point is the role of the boundary conditions exerted in the eastern Mediterranean: Experiments show that the initial geometry of the spreading layers partly controls their finite deformation and rotation and may explain the formation of several spreading arcs above the oceanic plate

Analogue Experiments

We use relatively simple analogue experiments simulating the Anatolian-Aegean tectonics, to study first-order tectonic phenomena in the area. For convenience, we shall refer to regions of the analogue experiment in terms of geographical directions (north, south, east, and west). The lithosphere in the Anatolian, Aegean, and European regions is modeled by a dry sand layer underlain by a silicone layer. The absence of cohesion and internal friction angle of the dry sand make it an analogue for the brittle upper crust, while the silicone putty simulates the ductile lower lithosphere (see the work of [START_REF] Davy | Experiments on shortening of 4-layer model of the continental lithosphere[END_REF] for details concerning scaling in this kind of experiments). These layers float above a liquid (honey) representing the asthenosphere, whose viscosity of about 20 Pa s, is negligible compared to the viscosity of silicone (Table 1). We did not include any details in the layers modeling the lithosphere to preserve the simplicity of the experiment. The sand and silicone layers are slightly lighter than the honey (1200 kg/m3 versus 1400 kg/m3; see Table 1) to simulate a floating continental lithosphere above the asthenosphere. Above the sand layer, a rectangular grid of passive colored-sand markers enables visualization of the deformation. The deformed grid is then digitized at different steps of the experiment to monitor the total and incremental displacements, deformation, and rigid rotations. The Mediterranean oceanic plate is simply modeled by honey without any other material above it. We are not interested here in examining the deformation of the oceanic plate but only in simulating the gravity spreading of the continental lithosphere above the subducting ocean [START_REF] Hatzfeld | Gravity spreading: A possible contribution to the geodynamics of the Aegean[END_REF].

Boundary Conditions

The sand-silicone plate is placed in a rectangular box (Figure 2). Tensional forces pull the southern boundary of the sand-silicone layer, because of the different buoyancy of the sand, silicone, and honey. When the forces arising from these density contrasts are large enough to break the sand layer, gravity spreading of the sand-silicone plate starts (Table 1). In nature the extension of the continental lithosphere above the retreating oceanic slab results from many other phenomena, as, for example, trench suction, which are not modeled here to preserve the simplicity of the experiments.

Crossing the box in the middle, a wall bounds the sand-silicone layer. At the beginning of the experiment, part of this wall is removed (Figure 2), and the light sand-silicone layer starts to spread over the honey, simulating the spreading of the continental lithosphere over the subduction zone. Note that the initial boundary of the spreading sand-silicone plate is not straight but angled at the eastern end.

Figure 2.

Boundary conditions for the analogue experiments. Thick lines represent rigid walls, and the thin line is the free boundary of the spreading sand-silicone layer.

A rigid indenter placed in the southeastern part of the box moves in the northward direction at a constant speed (a few centimenters per hour; see Table 1) to simulate the northward motion of Arabia with respect to stable Eurasia. On the other side of the box a small wall is placed to model the Apulian continental plate. In the experiment this wall does not move. We do not take into account the motion of Apulia with respect to Europe, because this motion is slower than the displacement of Arabia and also because Apulia is much smaller than Arabia. Thus it is assumed that the motion of Apulia is not responsible for the first-order structures that develop in the Anatolian-Aegean system. The other walls of the box are rigid. Their geometry does not have any influence on the deformation, since the deformed area is always confined in a small domain in the center of the box as [START_REF] England | Length scales for continental deformation in convergent, divergent, and strike-slip environments: Analytical and approximate solutions for a thin viscous sheet model[END_REF] predicted (see Figure 3). The advance of the rigid indentor starts at the beginning of the experiment, simultaneously with initiation of spreading of the sand-silicone layer. The speed of the piston is chosen to be roughly similar to the speed of the spreading sand-silicone layer over the honey, because in reality the speed of Arabia to the north is comparable with the speed of the Hellenic arc to the SW.

In section 4, we present only three different analogue experiments, although many more were built. The precise geometry of the faults that appear in the experiments we present, compared to similar ones we do not describe here, may not be exactly the same. The strain regime that develops in both experiments is, however, essentially the same, and the observations we discuss in section 5 have been reproduced in several experiments.

The Development of Lateral Escape in Experiments

As soon as the rigid wall constraining the sand-silicone layer is removed, spreading starts. Because of the initial geometry of the spreading layer, two arcs develop, both generated by several families of faults (Experiment 1, Figure 3). At the beginning of deformation during the first 20 min of Experiment 1, E-W normal faults develop near the E-W free boundary. Two families of conjugate strike-slip faults with a normal dip-slip component also develop within two oblique areas on both sides of this free boundary (Figure 3B). We observe a geometry of faults in this region similar to that obtained by [START_REF] Hatzfeld | Gravity spreading: A possible contribution to the geodynamics of the Aegean[END_REF] or [START_REF] Gautier | Timing, kinematics, and cause of Aegean extension: A scenario based on a comparison with simple analogue experiments[END_REF] in gravity spreading experiments without a moving piston. This result indicates that the westernmost arc that initially develops from the E-W boundary is not influenced by the advance of the piston. Between 20 and 70 min the shortened area in front of the piston grows northward. Dextral E-W strike-slip faults appear in the area situated between the shortened zone and the western arc (Figure 3). Figure 5 shows the incremental displacements of the nodes of the grid during this time period. Nodes situated in front of the piston move to the north, while nodes situated north of the eastern arc move to the west, and nodes situated within the western arc move to the southwest (Figure 4). This indicates that the shortening in front of the piston affects the whole deformed area and that extrusion starts during this time interval. The strike-slip faults which accommodate extrusion link preexisting normal faults of the stretched domain north of the western arc with preexisting thrusts situated in the shortened area in front of the piston (Figure 3). Strike-slip faults always appear after the thrust and normal faults. They result from the convergent growth of both normal and thrust faults. The position of strike-slip faults accommodating extrusion thus appears to be controlled by the location of the thrust and normal faults which developed in the shortened and stretched areas respectively.

To investigate this hypothesis, two different experiments were conducted. First, we conducted an experiment with the same sand and silicone densities as in Experiment 1 but with a larger free boundary (Experiment 2, Table 1). In this experiment the size of the stretched zone is very large compared to the size of the shortened area, and the orientation of dextral strike-slip faults is very different from that in Experiment 1, with lateral faults now oriented NW-SE rather than E-W, joining the northernmost grabens with the northernmost thrusts (Figure 6).

Second, we conducted Experiment 3, which has the same geometry as the reference experiment presented in Figures 3 and4. The density of the sand and silicone, however, is now close to the density of the honey substrate (Table 1). In this experiment, gravity spreading does not occur, because forces arising from lateral density contrasts are not strong enough to break the sand layer (Table 1). When the piston starts, shortening is concentrated in front of it without any lateral escape, even though a preexisting fault had been introduced in the sand layer to force extrusion (Figure 7). This experiment indicates that lateral escape cannot develop in the continental lithosphere if there is no extension somewhere to enable it [see also [START_REF] Suzanne | Extrusion latérale de l'Anatolie: géométrie et mécanisme de la fracturation[END_REF]. Using a finite element program to study the westward escape of Anatolia, [START_REF] Cianetti | Reproducing the velocity and stress fields in the Aegean region[END_REF] also concluded that extensional stresses in the Aegean region are necessary to enable the extrusion of Anatolia.

Figure 6. Experiment 2 after 145 min. The rheological layering of this experiment is the same

as that for Experiment 1 (Table 1). The main difference is the larger area where spreading occurs. Here the stretched zone is much more developed than the thrust area is. Dextral strikeslip faults accommodating lateral escape join the northernmost grabens with the thrust zone. These faults are oriented NW-SE in contrast to the E-W orientation observed in the experiment shown in Figure 3.

Figure 7. Experiment 3 after 76 min of shortening. The sand-silicone layer is heavier, and its density contrast with honey is low and does not exert enough force to break the sand layer and generate gravity spreading (Table 1). Shortening concentrates in front of the piston without any lateral escape, although a preexisting fault had been introduced in the model to facilitate extrusion.

As deformation induced by the piston movement continues (Experiment 1 after 83 min, Figure 3), the thrust zone in front of the piston propagates farther north of the stretched area.

The new compressional structures which appear at this time do not link anymore with the extensional grabens. Although the new thrusts faults accommodate a small strike-slip motion, lateral extrusion does not develop from them. There is no more lateral escape, and the spreading arc moves again to the south (Figures 3 and4).

Discussion

The Required Conditions for the Development of Extrusion

The most important implication of this experimental study is the indication that both the Arabia-Eurasia continental collision and the tensional deviatoric stresses arising from lateral density contrasts at the boundary between the Mediterranean oceanic plate and the Anatolian-Aegean continental domain, are likely requirements for the activation of the lateral escape of Anatolia (Figure 6) as was previously suggested by [START_REF] Suzanne | Extrusion latérale de l'Anatolie: géométrie et mécanisme de la fracturation[END_REF], [START_REF] Meijer | Temporal variation in the stress field of the Aegean region[END_REF], [START_REF] Cianetti | Reproducing the velocity and stress fields in the Aegean region[END_REF], and [START_REF] Lundgren | A geodynamic framework for eastern Mediterranean kinematics[END_REF]. The extrusion rate of Anatolia increases by approximately 10 mm/yr from east to west [START_REF] Reilinger | Global positioning system measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone[END_REF]. This increase indicates that the extrusion of the Anatolian plate probably is more controlled by the boundary conditions that prevail on the subduction zone than by the collision that occurs on the east. Experiments also indicate that the preexisting stretching of the Aegean lithosphere (which is older than the lateral escape of Anatolia [e.g., [START_REF] Gautier | Géométrie crustale et cinématique de l'extension tardi-orogénique dans le domaine centre-égéen; Ile des Cyclades et d'Eubée (Grèce)[END_REF]) may exert a strong control on the position of the dextral strike-slip faults which accommodate extrusion.

Deformation and Kinematics of Aegea and Anatolia

In Experiment 1 the kinematics and deformation associated with the lateral escape of the sand-silicone layer are very similar to the present-day kinematics of the Anatolian-Aegean region as deduced from geodetic data (compare Figures 1 and5) and seismicity. The western (Hellenic) arc moves to the SW as a whole, and within this arc the principal axes of the strain tensor are orientated roughly N-S (extension) and E-W (moderate shortening) (Figure 8). These directions correspond with those inferred from seismological or microtectonic data for the present-day deformation regime of the Aegean [e.g., [START_REF] Angelier | The tectonic development of the Hellenic arc and the Sea of Crete: a synthesis[END_REF][START_REF] Mercier | Extensional tectonic regimes in the Aegean basins during the Cenozoic[END_REF][START_REF] Hatzfeld | The strain pattern in the western Hellenic arc deduced from a microearthquake survey[END_REF][START_REF] Hatzfeld | Microearthquake seismicity and fault plane solutions in the southern Aegean and its tectonic implications[END_REF][START_REF] Hatzfeld | Microseismicity and strain pattern in northwestern Greece[END_REF].

Figure 8. Principal axes of the strain tensors calculated for triangles of the passive grid for Experiment 1 after 58 min of deformation. Thin lines indicate extension, and thick lines indicate shortening.

In the experiment the area situated north of the eastern arc (which corresponds in nature to central Anatolia) moves to the west. It does not show any N-S extension, because it corresponds to an intermediate domain situated between the collision zone imposed by the indentor (responsible for N-S shortening) and the western arc where gravity spreading prevails (responsible for N-S extension). The intermediate position of central Anatolia between Aegea and the Arabia-Eurasia collision zone could explain why N-S extension is not observed there (there is no N-S extension east of 31°E [START_REF] Angelier | The tectonic development of the Hellenic arc and the Sea of Crete: a synthesis[END_REF]). It may also explain why the Cypriot arc does not move to the south and why the activity of this arc is so moderate compared with the activity of the Hellenic arc although the Mediterranean oceanic lithosphere is present south of both the Aegean and Cypriot arcs. On one hand, Cyprus is pulled to the south by the oceanic lithosphere, while on the other, it is pushed to the north by the motion of the neighboring Arabian plate.

The Appearance of Two Distinct Spreading Arcs Above the Oceanic Plate.

Most authors [e.g., [START_REF] Meijer | Present-day dynamics of the Aegean region: A model analysis of the horizontal pattern of stress and deformation[END_REF] modeling the deformation of the Aegean region consider that the Aegean lithosphere spreads over the subducted Mediterranean oceanic lithosphere between the Keffalinia fault and the Antalya Gulf (Figure 1). The reason why the Aegean arc ends near Keffalinia is clear: The Keffalinia fault corresponds to the transition from the oceanic subduction to the continental collision with the Apulian plate [e.g., [START_REF] Kahle | The Kephalonia transform fault and the rotation of the Apulian platform: Evidence from satellite geodesy[END_REF][START_REF] Baker | Earthquake mechanisms of the Adriatic sea and western Greece: Implications for the oceanic subduction-continental collision transition[END_REF]. Apulia opposes the advance to the SW of the Aegean plate, and GPS data show that Epirus, north of the Keffalinia fault, is stable with respect to Europe [e.g., Le Pichon et al., 1995] (Figure 1). In contrast, the eastern edge of the Hellenic arc south of Antalya does not correspond to any rheological change of the Mediterranean plate. Seismic refraction data show that the Mediterranean oceanic plate is present to the east as far as south Cyprus and the South Levantine Basin [START_REF] Makris | Seismic refraction profiles between Cyprus and Israel and their interpretation[END_REF][START_REF] Makris | Physical properties and state of the crust and upper mantle of the Eastern Mediterranean Sea deduced from geological data[END_REF].

The boundary conditions that most workers adopt to reproduce the formation of the Antalya Gulf (no N-S motion of the southern boundary of the Anatolian plate east of Antalya; see, for instance, the work of [START_REF] Hatzfeld | Gravity spreading: A possible contribution to the geodynamics of the Aegean[END_REF] or [START_REF] Meijer | Present-day dynamics of the Aegean region: A model analysis of the horizontal pattern of stress and deformation[END_REF]) is convenient to get the eastern boundary of the Hellenic arc near Antalya. In reality, however, the same boundary conditions should apply south of the Aegean and Cypriot arcs. Then, such conditions cannot be invoked to explain the formation of the Antalya Gulf. Experiment 1 shows that the initial geometry of the continental plate, before the start of the subduction of the oceanic Mediterranean lithosphere, may have resulted in the formation of two distinct spreading arcs on both sides of the Antalya Gulf.

Rotations

The finite rigid body rotations evidenced in our experiments may be compared to paleomagnetic data. Clockwise rotations develop in the western parts of the western arc in Experiment 1. They reproduce the recent clockwise rotations measured using paleomagnetism in continental Greece [START_REF] Kissel | The Tertiary geodynamical evolution of the Aegean arc: A paleomagnetic reconstruction[END_REF] (Figure 1), although rotations are observed in the region north of the Keffalinia fault extending as far as Albania [START_REF] Speranza | Paleomagnetic evidence for Cenozoic clockwise rotation of the external Albanides[END_REF]. A possible explanation for the rotations measured north of the Keffalinia fault may be that oceanic or partially oceanic lithosphere was subducting till recently below Albania : The NW boundary of the Hellenic subduction would have migrated to the south, from northern Albania to Keffalinia. This may explain the recent deformation reported by [START_REF] Moretti | Deflection, gravity anomalies and tectonics of doubly subducted continental lithosphere: Adriatic and Ionian Seas[END_REF] in the Hellenides, south of the Scutari-Pec fault (northern Albania), and the apparent contradiction between paleomagnetic and geodetic data : The rapid motions responsible for the rotations reported in Albania by [START_REF] Speranza | Paleomagnetic evidence for Cenozoic clockwise rotation of the external Albanides[END_REF] for the Neogene (Mio-Pliocene) may have stopped with continental collision.

In Experiment 1, counterclockwise rotations appear in the eastern parts of the western arc, while no rotation occurs in the western part of the eastern arc. This is consistent with the absence of recent (Miocene) rotation in the Taurides east of the Antalya Gulf, while the western side of the gulf records a 30° counterclockwise rotation for the same period [START_REF] Kissel | The Tertiary geodynamical evolution of the Aegean arc: A paleomagnetic reconstruction[END_REF][START_REF] Kissel | First paleomagnetic evidence for a post-Eocene clockwise rotation of the western Taurides thrust belt east of the Isparta reentrant (southwestern Turkey)[END_REF]. It is clear, however, that the amount of rotation we observe in the experiment is highly dependent on the initial geometry adopted for the model. The pre-Miocene geometry of the Anatolian-Aegean plate may not be similar to the geometry that we adopted here. This may explain why the rotations we observe in Experiment 1 in the southern part of the eastern arc do not correspond with the recent rotations measured in Cyprus. Measurements by [START_REF] Abrahamsen | Palaeomagnetic timing of the rotation and translation of Cyprus[END_REF] indicate that Cyprus displaced 1500 km to the north during the last 30 Myr, but rotated counterclockwise less than 20° during the same period. Thus analogue experiments suggest that the different rotations measured on both sides of the Antalya Gulf may result from the initial geometry of the plate, which is responsible for the formation of the two spreading arcs.

The North Anatolian Fault

In experiments the dextral strike-slip motion which accommodates extrusion is generally distributed over several faults. In reality, however, the lateral escape of Anatolia is absorbed by the North Anatolian fault only. Le Pichon et al. [1995] show that the displacement vectors with respect to stable Europe for points of the Anatolian-Aegean plate are close to the vectors that would result from the rotation of a rigid plate around an Euler pole situated near the Nile delta. This means that the Anatolian-Aegean plate, as a first-order approximation, presently behaves as a rigid plate, namely, that its internal deformation is responsible for displacements smaller than the overall motion of the "plate" [START_REF] Reilinger | Global positioning system measurements of present-day crustal movements in the Arabia-Africa-Eurasia plate collision zone[END_REF]. The experiments presented here cannot reproduce this observation. We do not obtain in our experiment any "rigid" plate which behaves like Anatolia. The reason for this difference between nature and experiment may result from the heterogeneity of the rheology of the continental lithosphere. We know that major lithospheric faults such as the North Anatolian fault concentrate deformation since they are much weaker than the surrounding continental lithosphere [e.g., [START_REF] Tapponnier | On the mechanics of the collision between India and Asia[END_REF][START_REF] Rice | Fault stress states, pore pressure, distributions, and the weakness of the San Andreas Fault[END_REF]. In contrast, we have carried out experiments on homogeneous material to avoid introducing too many arbitrary assumptions. This suggests that the rigid plate behavior of Anatolia results from the weakness of the North Anatolian fault [START_REF] Cianetti | Reproducing the velocity and stress fields in the Aegean region[END_REF].

Conclusion

The experiments presented in this paper provide simple explanations for the present-day tectonics of the Anatolian-Aegean area, particularly for what concerns the existence of the two Hellenic and Cypriot arcs, and for the different rotations measured in paleomagnetism on each side of the Antalya Gulf [START_REF] Kissel | First paleomagnetic evidence for a post-Eocene clockwise rotation of the western Taurides thrust belt east of the Isparta reentrant (southwestern Turkey)[END_REF]. The experiments suggest that both the Arabia-Eurasia continental collision and the extension at the boundary between the European plate and the Mediterranean are necessary to generate the lateral extrusion of Anatolia in accordance with [START_REF] Suzanne | Extrusion latérale de l'Anatolie: géométrie et mécanisme de la fracturation[END_REF] or [START_REF] Cianetti | Reproducing the velocity and stress fields in the Aegean region[END_REF]. The geometry of the extruded plate may be controlled by the size of the thrust zone to the east and also by the size of the stretched area to the west.

Figure 3 .

 3 Figure 3. Top view of the reference experiment (Experiment 1; see Table 1) at different stages of deformation. (a) Beginning of experiment. The wall maintaining the sand-silicone layer has just been removed. (b) After 18 min of experiment. The spreading of the two arcs has started, and the shortening is concentrated very near the piston. Lateral escape still does not occur. (c) After 43 min of experiment. Dextral strike-slip faults accommodating the lateral extrusion have appeared. These faults join the thrusts in front of the piston with the grabens situated on the western side of the experiment. (d) After 83 min of experiment. The new thrusts in front of the piston appear too far north of the stretched zone to be connected to the grabens. There is no more lateral extrusion. (e) Line drawing of the main faults in the experiment after 43 min.

Figure 3 (

 3 Figure 3 (continued) : Figure 3e

Figure 4 .

 4 Figure 4. Displacement of passive markers on the western spreading arc in Experiment 1. At the beginning of the experiment (between 0 and 18 min), markers move to the south, showing that lateral escape has not started. During the interval 18-73 min, lateral escape prevails, and markers move to the southwest. At the end of the experiment, lateral escape has stopped, and markers move again to the south.

Figure 5 .

 5 Figure 5. Displacement of the nodes of the passive grid during the interval 33-58 min in Experiment 1.

Acknowledgments. Analogue experiments were built in the Modelisation Laboratory of Géosciences Rennes. Many thanks to marvelous J.J.K. Kermarrec for his technical and friendly assistance. C. Faccenna and J.P.J. Reijs provided careful and constructive reviews of the manuscript. P. Molnar and B. Chouet helped us to improve the manuscript.