Keywords: Modèle de Biot, ondes poroélastiques, solution analytique, technique de Cagniard de Hoop

Solution analytique pour la propagation d'ondes en milieu poroélastique stratifié. Partie I : en dimension 2

Résumé : Nous nous intéressons à la modélisation de la propagation d'ondes dans les milieux infinis bicouches poroélastiques. Nous considérons ici le modèle bi-phasique de Biot. Cette partie est consacrée au calcul de la solution analytique en dimension deux à l'aide de la technique de Cagniard-De Hoop. Les solutions que nous présentons ici ont été utilisés pour la validation de codes numériques.

Introduction

Many seismic materials cannot only be considered as solid materials. They are often porous media, i.e. media made of a solid fully saturated with a fluid: there are solid media perforated by a multitude of small holes (called pores) filled with a fluid. It is in particular often the case of the oil reservoirs. It is clear that the analysis of results by seismic methods of the exploration of such media must take to account the fact that a wave being propagated in such a medium meets a succession of phases solid and fluid: we speak about poroelastic media, and the more commonly used model is the Biot's model [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. higher frequency range[END_REF][START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF].

When the wavelength is large in comparison with the size of the pores, rather than regarding such a medium as an heterogeneous medium, it is legitimate to use, at least locally, the theory of homogeneization [START_REF] Burridge | Biot's poroelasticity equations by homogenization[END_REF][START_REF] Hornung | Homogenization and porous media[END_REF]. This leads to the Biot's model [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. higher frequency range[END_REF][START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF] which involves as unknown not only the displacement field in the solid but also the displacement field in the fluid. The principal characteristic of this model is that in addition to the classical P and S waves in a solid one observes a P "slow" wave, which we could also call a "fluid" wave: the denomination "slow wave" refers to the fact that in practical applications, it is slower (and probably much slower) than the other two waves.

The computation of analytical solutions for wave propagation in poroelastic media is of high importance for the validation of numerical computational codes or for a better understanding of the reflexion/transmission properties of the media. Cagniard-de Hoop method [START_REF] Cagniard | Reflection and refraction of progressive seismic waves[END_REF][START_REF] De Hoop | The surface line source problem[END_REF] is a useful tool to obtain such solutions and permits to compute each type of waves (P wave, S wave, head wave...) independently. Although it was originally dedicated to the solution to elastodynamic wave propagation, it can be applied to any transient wave propagation problem in stratified medium. However, as far as we know, few works have been dedicated to the application of this method to poroelastic medium. In [START_REF] Ezziani | Modélisation mathématique et numérique de la propagation d'ondes dans les milieux viscoélastiques et poroélastiques[END_REF] the analytical solution to poroelastic wave propagation in an homogeneous 2D medium is provided.

In order to validate computational codes of wave propagation in poroelastic media, we have implemented the codes Gar6more 2D [START_REF] Diaz | Gar6more 2d[END_REF] and Gar6more 3D [START_REF] Diaz | Gar6more 3d[END_REF] which provide the complete solution (reflected and transmitted waves) of the propagation of wave in stratified 2D or 3D media composed of acoustic/acoustic, acoustic/elastic, acoustic/poroelastic or poroelastic/poroelastic layers. The codes are freely downloadable at http://www.spice-rtn.org/library/software/Gar6more2D and http://www.spice-rtn.org/library/software/Gar6more3D.

We will focus in this paper on the 2D poroelastic case, the 2D acoustic/poroelastic case is detailed in [START_REF] Diaz | Analytical solution for wave propagation in stratified acoustic/porous media. part I: the 2D case[END_REF] and the three dimensional cases will be the object of forthcoming papers. The outline of the paper is as follows: we first present the model problem we want to solve and derive the Green problem from it (section 1). Then we present the analytical solution to the wave propagation problem in a two-layered 2D poroelastic (section 2). Finally we show how the analytical solution can be used to validate a numerical code (section 3).

The model problem

We consider an infinite two dimensional medium (Ω = R 2 ) composed of two homogeneous poroelastic layers Ω + = R×] -∞, 0] and Ω -= R × [0, +∞[ separated by an horizontal interface Γ (see Fig. 1). We first describe the equations in the two layers ( 1.1) and the transmission conditions on the interface Γ ( 1.2), then we present the Green problem from which we compute the analytical solution ( 1.3).

Ω + Ω -

First Layer

Second Layer

y = 0

Poroelastic equations

We consider the second-order formulation of the poroelastic equations [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range[END_REF][START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. II. higher frequency range[END_REF][START_REF] Biot | Mechanics of deformation and acoustic propagation in porous media[END_REF]:

                                   ρ Ü s + ρ f Ẅ -∇ ó Σ = F u , in Ω×]0, T ], ρ f Ü s + ρ w Ẅ + 1 K Ẇ + ∇P = F w , in Ω×]0, T ], Σ = λ∇ ó U s I 2 + 2µε(U s ) -β P I 2 , in Ω×]0, T ], 1 m P + β ∇ ó U s + ∇ ó W = F p , in Ω×]0, T ], U s (x, 0) = 0, W (x, 0) = 0, in Ω, U s (x, 0) = 0, Ẇ (x, 0) = 0, in Ω, (1)
with

(∇ ó Σ) i = 2 j=1 ∂Σ ij ∂x j ∀ i = 1, 2.
As usual I 2 is the identity matrix of M 2 (IR), and ε(U s ) is the solid strain tensor defined by:

ε ij (U ) = 1 2
∂U i ∂x j + ∂U j ∂x i .
In [START_REF] Biot | Theory of propagation of elastic waves in a fluid-saturated porous solid. I. low-frequency range[END_REF], the unknowns are:

U s the displacement field of solid particles;

W = φ(U f -U s )
, the relative displacement, U f being the displacement field of fluid particles and φ the porosity;

P , the fluid pressure;

Σ, the solid stress tensor.

The parameters describing the physical properties of the medium are given by:

ρ = φ ρ f + (1 -φ)ρ s
is the overall density of the saturated medium, with ρ s the density of the solid and ρ f the density of the fluid;

ρ w = aρ f /φ
, where a is the tortuosity of the solid matrix; K = κ/η, where κ is the permeability of the solid matrix and η is the viscosity of the fluid;

m and β are positive physical coefficients:

β = 1 -K b /K s and m = [φ/K f + (β -φ)/K s ] -1
, where K s is the bulk modulus of the solid, K f is the bulk modulus of the fluid and K b is the frame bulk modulus;

µ is the frame shear modulus, and λ = K b -2µ/3 is the Lamé constant.

F u , F w and F p are the force densities.

To simplify this study, we consider only the case of a compression source

F u (x, y, t) = f u ∇(δ x δ y-h ) f (t) and F w (x, y, t) = f w ∇(δ x δ y-h ) f (t)
and a pressure source F p = f p δ x δ y-h f (t), where f u , f w and f p are constant and f is a regular function source in time. We can generalize this approach for other types of punctual sources such as for instance

F u = f u ∇ × (δ x δ y-h ) f (t) and F w = f w ∇ × (δ x δ y-h ) f (t).

Transmission conditions

Let n be the unitary normal vector of Γ outwardly directed to Ω -. The transmission conditions on the interface Γ between the two poroelastic medium are [START_REF] Carcione | Wave Fields in Real Media : Wave propagation in Anisotropic, Anelastic and Porous Media[END_REF]:

           U + s = U - s , W + ó n = W -ó n, P + = P -, Σ + n = Σ -n.
(2)

The Green problem

We won't compute directly the solution to (2) but the solution to the following Green problem:

ρ ± ü± s + ρ ± f ẅ± -∇ ó σ ± = f u ∇(δ x δ y-h )δ t , in Ω ± ×]0, T ], (3a) ρ
± f ü± s + ρ ± w ẅ± + 1 K ± ẇ± + ∇p ± = f w ∇(δ x δ y-h )δ t , in Ω ± ×]0, T ], (3b)
σ ± = λ ± ∇ ó u ± s I 2 + 2µ ± ε(u ± s ) -β ± p ± I 2 , in Ω ± ×]0, T ], (3c)
1 m ± p ± + β ± ∇ ó u ± s + ∇ ó w ± = f p δ x δ y-h f (t), in Ω ± ×]0, T ], (3d)
u - s = u + , on Γ×]0, T ] (3e)
w -ó n = w + ó n, on Γ×]0, T ] (3f)

p -= p + , on Γ×]0, T ] (3g) σ -n = σ + n, on Γ×]0, T ]. (3h)
The solution to (1) is then computed from the solution to the Green Problem thanks to a convolution by the source function. For instance we have:

P + (x, y, t) = p + (x, y, .) * f (.) = t 0 p + (x, y, τ )f (t -τ ) dτ
(we have similar relations for the other unknowns). We also suppose that the poroelastic medium is non dissipative, i.e the viscosity η ± = 0. Using the equations (3c,3d) we can eliminate σ ± and p ± in (3) and we obtain the equivalent system:

         ρ ± ü± s + ρ ± f ẅ± -α ± ∇(∇ ó u ± s ) + µ ± ∇ × (∇ × u ± s ) -m ± β ± ∇(∇ ó w ± ) = (f u -β + m + f p )∇(δ x δ y-h ) δ t , ρ ± f ü± s + ρ ± w ẅ± -m ± β ± ∇(∇ ó u ± s ) -m ± ∇(∇ ó w ± ) = (f w -m + f p )∇(δ x δ y-h ) δ t , (4)
with

α -= λ -+ 2µ -+ m -β -2 .
And the transmission conditions on Γ are rewritten as:

                                         u + sx = u - sx , u + sy = u - sy , w - y = w + y , m + β + ∇ ó u + s + m + ∇ ó w + = m -β -∇ ó u - s + m -∇ ó w -, µ + (∂ y u + sx + ∂ x u + sy ) = µ -(∂ y u - sx + ∂ x u - sy ), (λ -+ m + β + 2 )∇ ó u + s + 2µ + ∂ y u + sy + m + β + ∇ ó w + = (λ -+ m -β -2 )∇ ó u - s + 2µ -∂ y u - sy + m -β -∇ ó w -. (5)
We split the displacement fields u ± s and w ± on irrotational and isovolumic fields (P-wave and S-wave):

u ± s = ∇Θ - u + ∇ × Ψ ± u ; w ± = ∇Θ ± w + ∇ × Ψ ± w . (6)
We can then rewrite system (4) in the following form:

                     A + Θ+ -B + ∆Θ + = δ x δ y-h δ t F , in Ω + ×]0, T ] A -Θ--B -∆Θ -= 0, in Ω -×]0, T ] Ψ± u -V ± S 2 ∆Ψ ± u = 0, in Ω ± ×]0, T ] Ψ± w = - ρ ± f ρ ± w Ψ± u , in Ω ± ×]0, T ] (7)
where

Θ ± = (Θ - u , Θ - w ) t , F = (f u -β + m + f p , f w -m + f p ) t , A ± and B ± are 2 × 2 symmetrical matrices: A ± = ρ ± ρ ± f ρ ± f ρ ± w ; B ± = λ ± + 2µ ± + m ± (β ± ) 2 m ± β ± m ± β ± m ± ,
and

V ± S = µρ ± w ρ ± ρ ± w -ρ ± f 2
is the S-wave velocity.

We multiply the first (resp. the second) equation of system [START_REF] De Hoop | The surface line source problem[END_REF] by the inverse of A + (resp. A -). The matrix

A + -1 B + (resp. A --1 B -) is diagonalizable: A ± -1 B ± = P ± D ± P ± -1
, where P ± is the change-of-coordinates matrix,

D ± = diag(V ± P f 2 , V ± P s
2 ) is the diagonal matrix similar to A ± -1 B ± , V ± P f and V ± P s are respectively the fast P-wave velocity and the slow Pwave velocity (V ± P s < V ± P f ).

Using the change of variables

Φ ± = (Φ ± P f , Φ ± P s ) t = P ± -1 Θ ± , (8)
we obtain the uncoupled system on fast P-waves, slow P-waves and S-waves:

                     Φ+ -D + ∆Φ + = δ x δ y-h δ t F + , in Ω + ×]0, T ] Φ--D -∆Φ -= 0, in Ω -×]0, T ] Ψ± u -V ± S 2 ∆Ψ ± u = 0, in Ω ± ×]0, T ] Ψ ± w = - ρ ± f ρ ± w Ψ ± u , in Ω ± ×]0, T ] (9)
with

F + = (A + P + ) -1 F = (F + P f , F + P s ) t .
Finally, we obtain the Green problem equivalent to (3):

                 Φ+ i -V + i 2 ∆Φ + i = δ x δ y-h δ t F + i , i ∈ {P f, P s} y > 0 Φ+ S -V + S 2 ∆Φ + S = 0 y > 0 Φ- i -V - i 2 ∆Φ - i = 0, i ∈ {P f, P s, S} y < 0 B(Φ + P f , Φ + P s , Φ + S , Φ - P f , Φ - P s , Φ - S ) = 0, y = 0 (10)
where we have set Φ ± S = Ψ ± u in order to have similar notations for the P f , P s and S waves. The operator B represents the transmission conditions on Γ:

B         Φ + P f Φ + P s Φ + S Φ - P f Φ - P s Φ - S         =                 P + 11 ∂ x P + 12 ∂ x ∂ y -P - 11 ∂ x -P - 12 ∂ x -∂ y P + 11 ∂ y P + 12 ∂ y -∂ x -P - 11 ∂ y -P - 12 ∂ y ∂ x P + 21 ∂ y P + 22 ∂ y ρ + f ρ + w ∂ x -P - 21 ∂ y -P - 22 ∂ y - ρ - f ρ - w ∂ x B 41 B 42 0 B 44 B 45 0 B 51 B 52 µ + (∂ 2 yy -∂ 2 xx ) B 54 B 55 -µ -(∂ 2 yy -∂ 2 xx ) B 61 B 62 -2µ + ∂ 2 xy B 64 B 65 2µ -∂ 2 xy                                Φ + P f Φ + P s Φ + S Φ - P f Φ - P s Φ - S               
INRIA where P ± ij , i, j = 1, 2 are the components of the change-of-coordinates matrix P ± and

B 41 = m + (β + P + 11 + P + 21 ) V + P f 2 ∂ 2 tt ; B 42 = m + (β + P + 12 + P + 22 ) V + P s 2 ∂ 2 tt ; B 44 = - m -(β -P - 11 + P - 21 ) V - P f 2 ∂ 2 tt ; B 45 = - m -(β -P - 12 + P - 22 ) V - P s 2 ∂ 2
B 64 = (λ -+ m -β -2 )P - 11 + m -β -P - 21 V - P f 2 ∂ 2 tt + 2µ -P - 11 ∂ 2 yy ; B 65 = (λ -+ m -β -2 )P - 12 + m -β -P - 22 V - P s 2 ∂ 2 tt + 2µ -P - 12 ∂ 2 yy .
To obtain this operator we have used the transmission conditions ( 5), the change of variables ( 6)-( 8) and the uncoupled system [START_REF] Diaz | Gar6more 2d[END_REF].

Moreover, from the unknowns Φ ± P f , Φ ± P s and Φ ± S we can determine the solid displacement u ± s and the relative displacement w ± by using the change of variables presented below.

Expression of the analytical solution

To state our results, we need the following notations and definitions:

1. Definition of the complex square root. For q ∈ C\IR -, we use the following definition of the square root g(q) = q 1/2 : g(q) 2 = q and ℜe[g(q)] > 0.

The branch cut of g(q) in the complex plane will thus be the half-line defined by {q ∈ IR -} (see Fig. 2). In the following, we use the abuse of notation g(q) = i √ -q for q ∈ IR -.

2. Definition of the functions κ ± i . For i ∈ {P f, P s, S} and q ∈ C, we define the functions

κ ± := κ + (q) = 1 V ± i 2 + q 2 1/2
.

ℑm(q)

ℜe(q) -π +π 3. Definition of the reflection and transmission coefficients. For a given q ∈ C, we denote by R P f P f (q), R P f P s (q), R P f S (q), T P f P f (q), T P f P s (q), and T P f S (q) the solution to the linear system

A(q)                 R P f P f (q)
R P f P s (q) R P f S (q)

T P f P f (q)

T P f P s (q)

T P f S (q)                 = 1 2κ + P f (q)V + P f 2                      i qP + 11 -κ + P f (q)P + 11 -κ + P f (q)P + 21 - m + V + P f 2 (β + P + 11 + P + 21 ) 2i qµ + κ + P f (q)P + 11 - (λ + + m + β + 2 )P + 11 + m + β + P + 21 V + P f 2 -2µ + P + 11 κ + P f (q)                     
and by R P sP f (q), R P sP s (q), R P sS (q), T P sP f (q), T P sP s (q) and T P sS (q) the solution to the linear system

A(q)                 R P sP f (q) R P sP s (q)
R P sS (q)

T P sP f (q) T P sP s (q)

T P sS (q)                 = 1 2κ + P s (q)V + P s 2                     i qP + 12 -κ + P s (q)P + 12 -κ + P s (q)P + 22 - m + V + P s 2 (β + P + 12 + P + 22 ) 2i qµ + κ + P s (q)P + 12 - (λ + + m + β + 2 )P + 12 + m + β + P + 22 V + P s 2 -2µ + P + 12 κ + P s (q)                    
, where the matrix A(q) is defined for q ∈ C by:

A(q) =                  -i qP + 11 -i qP + 12 -κ + S (q) i qP - 11 i qP - 12 -κ - S (q)
-κ + P f (q)P + 11 -κ + P s (q)P + 12 i q -κ - P f (q)P - 11 -κ - P s (q)P -

12

-i q -κ + P f (q)P + 21 -κ + P s (q)P +

22

-i q

ρ + f ρ + w -κ - P f (q)P - 21 -κ - P s (q)P - 22 i q ρ - f ρ - w A 41 (q)
A 42 (q) 0 A 44 (q) A 45 (q) 0 A 51 (q) A 52 (q) A 53 (q) A 54 (q) A 55 (q) A 56 (q)

A 61 (q) A 62 (q) -2i qµ + κ + S (q) A 64 (q)

A 65 (q) -2i qµ -κ - S (q)                  , with A 41 (q) = m + V + P f 2 β + P + 11 + P + 21 ; A 42 (q) = m + V + P s 2 β + P + 12 + P + 22 ; A 44 (q) = - m - V - P f 2 β -P - 11 + P - 21 ; A 45 (q) = - m - V - P s 2 β -P - 12 + P - 22 ;
A 51 (q) = 2i qµ + κ + P f (q)P + 11 ; A 52 (q) = 2i qµ + κ + P s (q)P + 12 ; A 53 (q) = µ + (κ + S 2 (q) + q 2 );

A 54 (q) = 2i qµ -κ - P f (q)P - 11 ; A 55 (q) = 2i qµ -κ - P s (q)P - 12 ; A 56 (q) = -µ -(κ - S 2 (q) + q 2 );

A 61 (q) = (λ + + m + β + 2 )P + 11 + m + β + P + 21

V + P f 2 + 2µ + κ + P f
2 (q)P + 11 ;

A 62 (q) = (λ + + m + β + 2 )P + 12 + m + β + P + 22

V + P s 2 + 2µ + κ + P s (q) 2 P + 12 ; A 64 (q) = - (λ -+ m -β -2 )P - 11 + m -β -P - 21 V - P f 2 -2µ -κ - P f 2 (q)P - 11 ; A 65 (q) = - (λ -+ m -β -2 )P - 12 -m -β -P - 22 V - P s 2 + 2µ -κ - P s 2 (q)P - 12 .
We also denote by V max the greatest velocity in the medium:

V max = max(V + P f , V + P s , V + S , V - P f , V - P s , V - S ).
We can now present the expression of the solution to the Green Problem:

Theorem 2.1. The solid displacement in the top medium is given by

u + s (x, y, t) = t 0 ν + (x,
             ν + P f,x (x, y, t) = - P + 11 F + P f V + P f 2 tx 2πr t 2 -t 2 0 , ν + P f,y (x, y, t) = - P + 11 F + P f V + P f 2 t(y -h) 2πr t 2 -t 2 0 , if t > t 0 ν + P f (x, y, t) = 0. else
We set here r = (x 2 + (y -h) 2 ) 1/2 and t 0 = r/V + P f denotes the arrival time of the incident P f wave. ν + P s is the velocity of the incident P s wave and satisfies:

           ν + P s,x (x,
V + P s 2 t(y -h) 2πr t 2 -t 2 0 , if t > t 0 ν + P s (x, y, t) = 0. else
We set here r = (x 2 + (y -h) 2 ) 1/2 and t 0 = r/V + P s denotes the arrival time of the incident P s wave.

ν +

P f P f is the velocity of the reflected P f P f wave (the P f reflected wave generated by the P f incident wave) and satisfies:

               ν + P f P f,x (x, y, t) = ℑm i υ(t)κ + P f (υ(t))R P f P f (υ(t)) π t 2 0 -t 2 P + 11 F + P f , ν + P f P f,y (x, y, t) = ℑm κ + P f 2 (υ(t))R P f P f (υ(t)) π t 2 0 -t 2 P + 11 F + P f , if t h < t ≤ t 0 and x r > V + P f V max INRIA                ν + P f P f,x (x, y, t) = - ℜe i γ(t)κ + P f (γ(t))R P f P f (γ(t)) π t 2 0 -t 2 P + 11 F + P f , ν + P f P f,y (x, y, t) = - ℜe κ + P f 2 (γ(t))R P f P f (γ(t)) π t 2 -t 2 0 P + 11 F + P f , if t > t 0 ν + P f (x, y, t) = 0. else
We set here r = (x 2 + (y + h) 2 ) 1/2 and t 0 = r/V + P f denotes the arrival time of the reflected P f P f volume wave,

t h = (y + h) 1
V + P f 2 - 1 V 2 max + |x| V max
denotes the arrival time of the reflected P f P f head wave and the complex functions υ := υ(t) and γ := γ(t) are defined by

                     υ(t) = -i   y + h r 1 V + P f 2 - t 2 r 2 + xt r 2   for t h < t ≤ t 0 and x < 0, υ(t) = i   y + h r 1 V + P f 2 - t 2 r 2 - xt r 2   for t h < t ≤ t 0 and x ≥ 0, and
γ(t) = -i x r 2 t + y + h r t 2 r 2 - 1 V + P f 2
for t > t 0 .

ν + P f P s is the velocity of the reflected P f P s wave and satisfies: ν + P f P s (x, y, t) = 0. else Here t 0 denotes the arrival time of the reflected P f P s wave (its calculation is similar to the calculation of the arrival time of the transmitted wave, see the appendix of [START_REF] Diaz | Analytical solution for wave propagation in stratified acoustic/porous media. part I: the 2D case[END_REF]) and t h denotes the arrival time of the P f P s head wave:

             ν + P f P s,x (x, y, t) = - P + 12 F + P f π ℜe i υ(t)R P f P s (υ(t)) dυ dt (t) , ν + P f P s,y (x, y, t) = - P + 12 F + P f π ℜe κ + P s (υ(t))R P f P s (υ(t)) dυ dt (t) , if t h < t ≤ t 0 and |ℑm (υ(t 0 ))| > 1 V max            ν + P f P s,x (x,
t h = y 1 V + P s 2 - 1 V 2 max + h 1 V + P f 2 - 1 V 2 max + |x| V max .
For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of q ∈ C → F(q, t) = y 1

V + P s 2 + q 2 1/2 + h 1
V + P f 2 + q 2 1/2 + iqx -t,
such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive. ν + P f S is the velocity of the reflected P f S wave and satisfies:

             ν + P f S,x (x, y, t) = - F + P f π ℜe κ + S (υ(t))R P f S (υ(t)) dυ dt (t) , ν + P f S,y (x, y, t) = F
+ P f π ℜe i υ(t)R P f S (υ(t)) dυ dt (t) , if t h < t ≤ t 0 and |ℑm (υ(t 0 ))| > 1 V max            ν + P f S,x (x, y, t) = - F + P f π ℜe κ + S (γ(t))R P f S (γ(t)) dγ dt (t) , ν + P f S,y (x, y, t) = F + P f π ℜe i γ(t)R P f S (γ(t)) dγ dt (t) , if t > t 0 ν + P f S (x, y, t) = 0. else
Here t 0 denotes the arrival time of the reflected P f S wave and t h denotes the arrival time of the reflected P f S head wave:

t h = y 1 V + S 2 - 1 V 2 max + h 1 V + P f 2 - 1 V 2 max + |x| V max .
For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

q ∈ C → F(q, t) = y 1 V + S 2 + q 2 1/2 + h 1 V + P f 2 + q 2 1/2 + iqx -t,
such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive. ν + P sP f is the velocity of the reflected P sP f wave and satisfies:

           ν + P sP f,x (x, y, t) = - P + 11 F + P s π ℜe i υ(t)R P sP f (υ(t)) dυ dt (t) , ν + P sP f,y (x, y, t) = - P + 11 F + P s π ℜe κ + P f (υ(t))R P sP f (υ(t)) dυ dt (t) , if t h < t ≤ t 0 and |ℑm (υ(t 0 ))| > 1 V max          ν + P sP f,x (x, y, t) = - P + 11 F + P s π ℜe i γ(t)R P sP f (γ(t)) dγ dt (t) , ν + P sP f,y (x, y, t) = - P + 11 F + P s π ℜe κ + P f (γ(t))R P sP f (γ(t)) dγ dt (t) , if t > t 0 ν + P sP f (x, y, t) = 0. else
Here t 0 denotes the arrival time of the reflected P sP f wave and t h denotes the arrival time of the reflected P sP f head wave:

t h = y 1 V + P f 2 - 1 V 2 max + h 1 V + P s 2 - 1 V 2 max + |x| V max .
For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of q ∈ C → F(q, t) = y 1

V + P f 2 + q 2 1/2
+ h 1

V + P s 2 + q 2 1/2 + iqx -t,
such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive. ν + P sP s is the velocity of the reflected P sP s wave and satisfies: We set here r = (x 2 + (y + h) 2 ) 1/2 and t 0 = r/V + P s denotes the arrival time of the reflected P sP s volume wave, t h = (y + h) 1

               ν + P sP s,x (x,
V + P f 2 - 1 V 2 max + |x| V max
denotes the arrival time of the reflected P sP s head wave and the complex functions υ := υ(t) and γ := γ(t) are defined by

             υ(t) = -i y + h r 1 V + P s 2 - t 2 r 2 + xt r 2
for t h < t ≤ t 0 and x < 0,

υ(t) = i y + h r 1 V + P s 2 - t 2 r 2 -
t h = y 1 V + S 2 - 1 V 2 max + h 1 V + P s 2 - 1 V 2 max + |x| V max .
For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

q ∈ C → F(q, t) = y 1 V + S 2 + q 2 1/2 + h 1 V + P s 2 + q 2 1/2 + iqx -t,
such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

Remark 2.1. For the practical computations of the displacement, we won't have to explicitly compute the primitive of the velocities ν, which would be rather tedious, since

t 0 ν(τ ) dτ * f = ν * t 0 f (τ ) dτ .
Therefore, we'll only have to compute the primitive of the source function f . The proof of this theorem is similar to the one detailed in [START_REF] Diaz | Analytical solution for wave propagation in stratified acoustic/porous media. part I: the 2D case[END_REF] for the computation of the analytical solution to the acoustic/poroelastic problem, therefore we won't detail it here.

Numerical illustration

To illustrate the use our results, we have compared our analytical solution to a numerical one obtained by C. Morency and J. Tromp [START_REF] Morency | Spectral-element simulations ofwave propagation in porous media[END_REF]. We consider an two-layered poroelastic medium whose characteristic coefficients are the solid density: The source is located in the acoustic layer, at 500 m from the interface. We used two types of sources in space: the first one is a bulk source such that f u = f w = -10 10 and f p = 0; the second one is a pressure source such that f u = f w = 0 and f p = 1. In each case we used a fifth derivative of a Gaussian of dominant frequency f 0 = 15 Hz:

f (t) = 4 π 2 f 2 0 9 t - 1 f 0 + 4 π 2 f 2 0 t - 1 f 0 3 -4 π 4 f 4 0 t - 1 f 0 5 e -π 2 f 2 0 t-1 f 0 2
for the source in time. We compute the solution at two receivers, the first one is in the upper layer, at 533 m from the interface; the second one is in the bottom layer, at 533 m from the interface; both are located on a vertical line at 400 m from the source (see Fig. 3). We represent the y component of the displacement from t = 0 to t = 1 s in Fig. 4 for the bulk source and in Fig. 5 for the pressure source. The left pictures represents the solution at receiver 1 while the right pictures represents the solution at receiver 2. On all the pictures the blue solid curve is the analytical solution and the red dashed curve is the numerical solution.

All the pictures show a good agreement between the two solutions, which validates the numerical code.

Conclusion

We provided the complete solution (reflected and transmitted wave) of the propagation of wave in a two-layered 2D poroelastic medium and we used it to validate a numerical code.

In a forthcoming paper we will use this solution as a basis to derive the solution in a three dimensional medium.

Figure 1 :

 1 Figure 1: Configuration of the study

Figure 2 :

 2 Figure 2: Definition of the function x → (x) 1/2

  ρ + s = 2200 kg/m 3 and ρ - s = 2650 kg/m 3 ; the fluid density: ρ + f = 950 kg/m 3 and ρ - f = 750 kg/m 3 ;the porosity: φ + = 0.4 and φ -= 0.2 ;the tortuosity: a + = 2 and a -= 2;the solid bulk modulus: K + s = 6.9 GPa and K - s = 37 GPa;the fluid bulk modulus:K + f = 2GPa and K - f = 1.7 GPa; the frame bulk modulus: K + b = 6.7 GPa and K - b = 2.2 GPa; the frame shear modulus µ + = 3 GPa and µ -= 4.4 GPa; so that the celerity of the waves in the poroelastic medium are: for the fast P wave, V + P f = 2692 m/s and V - P f = 2535 m/s; for the slow P wave, V + P s = 1186 m/s and V - P s = 744 m/s; for the ψ wave, V + S = 1409 m/s and V - S = 1415 m/s.

Figure 4 :

 4 Figure 3: Configuration of the experiment

Figure 5 :

 5 Figure5: The y component of the displacement at receiver 1 (left picture) and 2 (right picture) in the case of a pressure source. The blue solid curve is the analytical solution computed by the Cagniard-de Hoop method, the red dashed curve is the numerical solution.

  denotes the arrival time of the reflected P sS wave and t h denotes the arrival time of the reflected P sS head wave:

								xt r 2	for t h < t ≤ t 0 and x ≥ 0,
	and
			γ(t) = -i	x r 2 t +	y + h r	t 2 r 2 -	1 P s V +	2	for t > t 0 .
	ν + P sS is the velocity of the reflected P sS wave and satisfies:
	     	ν + P sS,x (x, y, t) = -	F + P s π	ℜe κ + S (υ(t))R P sS (υ(t))	dυ dt	(t) ,	if t h < t ≤ t 0
	    	ν + P sS,y (x, y, t) =	F + P s π	ℜe i υ(t)R P sS (υ(t))	dυ dt	(t) ,	and |ℑm (υ(t 0 ))| >	1 V max
	        	ν + P sS,x (x, y, t) = -ν + P sS,y (x, y, t) = F + F + P s π P s π ℜe i γ(t)R P sS (γ(t)) ℜe κ + S (γ(t))R P sS (γ(t)) dt (t) , dγ dt dγ	(t) ,	if t > t 0
	ν + P sS (x, y, t) = 0.							else
	Here t 0

RR n 6591

Acknowledgments

We thanks Christina Morency who provided us the numerical solutions we have used to validate our analytical solution.

ν - P f P f is the velocity of the transmitted P f P f wave (the P f transmitted wave generated by the P f incident wave) and satisfies:

Here t 0 denotes the arrival time of the transmitted P f P f wave and t h denotes the arrival time of the transmitted P f P f head wave:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of q ∈ C → F(q, t) = -y 1

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

ν - P f P s is the velocity of the transmitted P f P s wave and satisfies: ν - P f P s (x, y, t) = 0. else Here t 0 denotes the arrival time of the transmitted P f P s wave and t h denotes the arrival time of the transmitted P f P s head wave:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

ν - P f S is the velocity of the transmitted P f S wave and satisfies:

Here t 0 denotes the arrival time of the transmitted P f S wave and t h denotes the arrival time of the transmitted P f S head wave:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

ν -

P sP f is the velocity of the transmitted P sP f wave and satisfies:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of q ∈ C → F(q, t) = -y 1

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

ν -

P sP s is the velocity of the transmitted P sP s wave and satisfies:

ν - P sP s (x, y, t) = 0. else Here t 0 denotes the arrival time of the transmitted P sP s wave and t h denotes the arrival time of the transmitted P sP s head wave:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

+ h 1

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.

ν - P sS is the velocity of the transmitted P sS wave and satisfies:

For t h < t ≤ t 0 , the function υ(t) is implicitly defined as the only root of

+ h 1

such that ℑm x dυ(t) dt ≤ 0. For t > t 0 , the function γ(t) is defined as the only root of q ∈ C → F(q, t) whose real part is positive.


Bruno Martin

Patrick Solé

Pseudo-random Sequences Generated by Cellular Automata

Generation of pseudo random sequences by cellular automata, as well as by hybrid cellular automata is surveyed. An application to the fast evaluation and FPGA implementation of some classes of boolean functions is sketched out.

Notations and definitions 1.(Pseudo-)randomness

This section recalls the classical definitions of pseudo-randomness. We first give an intuitive statement which gives the difference between real randomness and pseudorandomness. We then introduce more formal definitions of pseudo-randomness.

Introduction

Cellular Automata (CA) is a popular model of finite state machine with some pretention to generality and universality. Pseudo Random Sequences (PRS) on the other hand, have a long history of applications to computational (Monte Carlo sampling, numerical simulation) and comunications problems (coding theory, streamciphers). In that context the popular model is the Linear Feedback Shift Register (LFSR), another model of linear finite state machine.

In the present work we survey the known attempts to generate PRS by CA. We give an account of the synthesis of LFSR by arrays of variable CA (known as hybrid CA or HCA). We sketch an application to the evaluation of boolean functions in n variables which are related to cyclic codes of length 2 n -1. This is aimed at VLSI implementation, especially by programmable arrays.

The material is organized as follows. Section 1 collects definitions and basic notions on PRS, CA and HCA. Section 2 reviews the synthesis theory of LFSR by HCA. Section 3 surveys the generation of PRS by elementary CA. Section 4 surveys the generation of PRS by HCA. Section 5 contains the application of synthesis theory to boolean functions evaluation.

In [START_REF] Wolfram | A new kind of science[END_REF], Wolfram describes three mechanisms responsible for random behavior in systems: [START_REF] Augot | On the computation of minimal polynomials, cyclic vectors, and frobenius forms[END_REF] Randomness from physics like brownian motion; (2) Randomness from the initial conditions which is studied by chaos theory; and (3) Randomness by design, also called pseudo-randomness used in pseudo-random sequences generators. Many algorithms generate pseudo-random sequences. The behavior of the system is fully determined by knowing the seed and the algorithm used. They are quicker methods than getting "true" randomness from the environment, inaccessible for computers.

The applications of randomness have led to many different methods for generating random data. These methods may vary as to how unpredictable or statistically random they are, and how quickly they can generate random sequences. Before the advent of computational PRS, generating large amount of sufficiently random numbers (important in statistics and physical experimentation) required a lot of work. Results would sometimes be collected and distributed as random number tables or even CD iso-images.

More formally, a pseudo-random sequence (PRS for short) can be defined as: Definition 1. A sequence is pseudo-random if it cannot be distinguished from a truly random sequence by any efficient (polynomial time) procedure or circuit.

Theorem 1 ([2]

). A sequence is pseudo-random iff it is next-bit unpredictable.

Theorem 1 claims that for pseudo-random sequences, even if we know all the history, we don't have any information on the next bit. Theorem 1 was proved equivalent to:

Theorem 2 ([22]).

A PRS generator G passes Yao's test if, for any family of circuits F with a polynomial number of gates for computing a statistical test, G passes F .

Cellular automata

In this section, we recall several definitions of cellular automata (CA). We focus on elementary cellular automata rules which restrict the set of the states to be F 2 . A cellular automaton is generally a bi-infinite array of identical cells which evolve synchronously and in parallel according to a local transition function. The cells can only communicate with their nearest neighbors. Here, we will concentrate on two finite restrictions of CA:

-cyclic: a ring of N cells indexed by Z N .

-null boundary: an array of N cells in which both extremal cells are fed with zeroes.

All the cells are finite state machines with a finite number of states and a transition function which gives the new state of a cell according to its current state and the current states of its nearest neighbors.

Definition 2. A cellular automaton is a finite array of cells. Each cell is a finite state machine

C = (Q, f ) where Q is a finite set of states and f a mapping f : Q 3 → Q.
The mapping f , called local transition function, has the following meaning: the state of cell i at time t + 1 (denoted by x t+1 i ) depends upon the state of cells i -1, i and i + 1 at time t (the neighborhood of cell i of radius 1). Fig. 1 illustrates one transition of a cellular automaton with 8 cells. The following equality rules the dynamics of the cellular automaton:

x t+1 i = f (x t i-1 , x t i , x t i+1 ) (1)
0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 t=0 t=1 neighborhood 0 1 7
Fig. 1. Transition of a cell (rule 30); cyclic CA.

For a fixed t, the sequence of all the values x t i for i ∈ Z N , is the configuration at time t. It is a mapping c which assigns a state of Q to each cell of the cellular automaton. The sequence of configurations as pictured by Fig. 2 is called a time-space diagram. Fig. 2 depicts the evolution of a ring with N = 8 cells. On the top of Fig. 2, we have depicted rule 30 with each transition illustrated by three adjacent squares representing the different preimages of f and on the bottom, their image by f . A 0 (resp. 1) is painted white (resp. black). On the bottom of Fig. 2, we see the time-space diagram of the cellular automaton from the initial configuration at time t = 0 to time t = 7.

LFSR synthesis by HCA

We will restrict ourselves to the case where Q = F 2 and f is a Boolean predicate with 3 variables, an elementary rule. These CA have been considered in [START_REF] Wolfram | Theory and applications of cellular automata[END_REF]: there are 256 different binary CA and a natural number can be associated to each rule as follows:

x t i-1 x t i x t i+1
111 110 101 100 011 010 001 000

x t+1 i 0 0 0 1 1 1 1 0
The top line gives all possible preimages for f and the bottom line the images by f . Thus, f is fully specified by the 8-bit number written on the bottom line (00011110 in our example) which can be translated in basis 10 and then called the rule of the cellular automaton (as rule number 30 here). Equivalently, this rule can be considered as a Boolean function with (at most) 3 variables. Taking rule 30 again, its corresponding Boolean function is:

x t+1 i = x t i-1 ⊕ (x t i ∨ x t i+1 )
with ⊕ denoting the Boolean XOR function and ∨ the classical Boolean OR function. Its equivalent formulation in F 2 is: Equivalent rules Since we are dealing with pseudo-random generators, some of the elementary rules are equivalent by three transformations, all introduced by Wolfram in [19, p. 492]. We first introduce some notation: let us denote by w the mirror image of the finite binary word w = w 1 . . . w n , w = w n . . . w 1 and by w the word obtained from w by exchanging the 0's by 1's (and conversely) w = w 1 . . . w n . The first transformation is the conjugation which interchanges the roles of 0 and 1. It takes as an input r, the binary representation of a rule and returns r. For instance, the conjugation transforms rule 30 into rule 135. The second transformation, called reflection gives a re-ordering of the bits of r.

x t+1 i = x t i-1 + x t i + x t i+1 + x t i x t i+1 . t=0
Each bit f r (x i-1 , x i , x i+1 ) is replaced by the value of f r (x i+1 , x i , x i-1 ) (the mirror image of x i-1 x i x i+1
) and leads to a re-ordering of the bits of r, the binary representation of the rule. As an example, by reflection, rule 30 is changed into rule 86. The last transformation combines boths and is called conjugation-reflection; it changes rule 30 into rule 149. All of these transformations keep the Walsh-Hadamard transform values of the cellular automata dynamics and are thus statistically equivalent.

Hybrid Cellular Automata

In the sequel, we will consider the case where different cells of the CA can use different rules. This model is called hybrid and will be denoted by HCA for short. In the context of sequence generation, several authors have considered this extension of the model of CA [START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF][START_REF] Seredynski | Cellular automata computations and secret key cryptography[END_REF][START_REF] Sipper | Co-evolving parallel random number generators[END_REF]]. We will focus on linear HCA (LHCA) which is used by [START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF].

Linear hybrid cellular automata In [START_REF] Serra | The analysis of one dimensional CA and their aliasing properties[END_REF][START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF]

M = [d 0 , d 1 , . . . , d N -1 ] such that d i = 0 if cell i uses rule 90 1 if cell i uses rule 150 . Given M , its reversal is M 's mirror image: [d N -1 , . . . , d 1 , d 0 ]
. We also define the subvector M i,j = [d i , . . . , d j ] with i ≤ j which also represents a submachine of the HCA consisting of cells i through j.

The encoding of rules 90 and 150 into zero and one, resp., means that equation ( 1) can be rewritten in F 2 as

x t+1 i = f i (x t i-1 , x t i , x t i+1 ) = x t i-1 + d i x t i + x t i+1
. We define the state of a HCA at time t to be the n-tuple formed from the state of the cells:

x t = [x t 0 , x t 1 , . . . , x t N -1 ]
T (the superscript T denotes the transpose). Then, the next state function of the HCA is computed as x t+1 = f (x t ). Since each f i is linear, f is also linear and an endomorphism of F N 2 . Linearity implies the existence of a matrix A such that x t+1 = f (x t ) = A • x t . The HCA transition matrix plays the same role as an LFSR transition matrix. A is tridiagonal.

A =        d0 1 0 ••• ••• 0 0 1 d1 1 . . . 0 0 1 d2 . . . . . . . . . . . . 1 dN-2 1 0 0 ••• ••• 0 1 dN-1        Let us denote by ∆ the characteristic polynomial of A, that is ∆ =| xId -A |.

Definition 3. [4]

A polynomial p is said to be a HCA polynomial if it is the characteristic polynomial of some HCA.

Recall that M i,j is the HCA consisting of cells i through j and denote ∆ i,j its corresponding characteristic polynomial. When i = 0, we simply write M k (resp. ∆ k ) for the CA consisting of cells 0 to k (resp. its corresponding characteristic polynomial). Cattell and Muzio [START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF] proved that ∆ k satisfies a recurrence relation:

Theorem 3. [4] ∆ k satisfies the reccurrence: ∆ -2 = 0, ∆ -1 = 1, ∆ k = (x + d k )∆ k-1 + ∆ k-2 for k ≥ 0.
Theorem 3 provides an efficient algorithm to compute ∆ N -1 the characteristic polynomial of a HCA from its rule vector M . Actually, this recurrence relation is related to Euclidean GCD algorithm on polynomials with ∆ k as the dividend, ∆ k-1 as the divisor, x + d k as the quotient and ∆ k-2 as the remainder. Applying Euclid's extended greatest division algorithm yields to the sequence of quotients whose constant terms are the mirror image of the rule vector. This comes from:

Lemma 1. [4] Let p ∈ F 2 [x] and q ∈ F 2 [x]
of respective degrees n and n -1. Then there exists a HCA with characteristic polynomial p and characteristic subpolynomial q if and only if applying Euclid's greatest division algorithm to p and q results in n degree one quotients.

Thus ∆ N -1 and ∆ N -2 determine the whole HCA. But in general, a characteristic polynomial isn't sufficient to uniquely determine the HCA. Just consider the following counter-example: [0, 0, 1, 0, 0, 0] ↔

x 6 + x 5 + x 4 + x 3 + 1 ↔ [1, 1, 0, 1, 1, 1]
To uniquely determine the HCA, we must know one more characteristic subpolynomial ∆ 1,N -1 and use theorem 4: Theorem 4 (HCA quadratic congruence [START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF]). Suppose we have a HCA with characteristic polynomial ∆ N -1 and characteristic subpolynomials ∆ N -2 and ∆ 1,N -1 . Then both y = ∆ N -2 and y = ∆ 1,N -1 satisfy the congruence:

y 2 +(x 2 +x)∆ ′ N -1 y +1 ≡ 0 mod ∆ N -1 where ∆ ′ N -1 is the formal derivative of ∆ N -1 in F 2 .
By combining Lemma 1 and Theorem 4, Cattell and Muzio give a characterization of HCA polynomials and give an algorithm for finding a HCA given a polynomial. Their method has been recently improved in [START_REF] Cho | New synthesis of onedimensional 90/150 linear hybrid group CA[END_REF].

Corollary 1. Let p ∈ F 2 [x] of degree n.
Then p is a HCA polynomial if and only if for some solution q for y of the congruence

y 2 + (x 2 + x)p ′ y + 1 ≡ 0 mod p (2)
Euclid's greatest division algorithm on p and q results in n degree one quotients.

Theorem 4 has some weaknesses: it does not say neither that polynomials solutions to the quadratic congruence will be subpolynomials of ∆ N -1 nor that non HCA polynomials won't have solutions to the quadratic congruence. Theorem 4 only gives necessary conditions for HCA polynomials: they have solutions to the quadratic congruence and that some of these solutions are subpolynomials. However, Theorem 4 is useful for irreducible polynomials:

Theorem 5. If p ∈ F 2 [x]
is an irreducible polynomial of degree n, then equation ( 2) has exactly two solutions, both of which result in n degree one quotients.

Corollary 2. If p ∈ F 2 [x]
is an irreducible polynomial, then p has exactly two HCA realizations with one being the reversal of the other.

Since one can build a HCA from an irreducible polynomial and represent it by its transition matrix, we can ask which is the relationship between LHCA and LFSR. If both are based on the same irreducible or primitive polynomial, they have the same behavior up to permutation of the order in which the states appear and the cycle structure of the states is identical. A similarity transform between LHCA and LFSR has been given in [START_REF] Cattell | An explicit similarity transform between CA and LFSR matrices[END_REF] and recently improved in [START_REF] Kagaris | A similarity transform for linear finite state machines[END_REF].

In [START_REF] Wolfram | Cryptography with CA[END_REF][START_REF] Wolfram | Random sequence generation by cellular automata[END_REF], Wolfram uses a one-dimensional cellular automaton for pseudo-random bit generation by selecting the values taken by a single cell when iterating the computation of rule 30 from an initial finite configuration where the cells are arranged on a ring of N cells. Mathematically, Wolfram claims the sequence {x t i } t≥0 is pseudo-random for a given i. Wolfram extensively studied this particular rule, demonstrating its suitability as a high performance randomizer which can be efficiently implemented in parallel; indeed, this is one of the pseudo-random generators which was shipped with the connection machine CM2 and which is currently used in the Mathematica R software.

Unfortunately, this PRG is not suitable for cryptographic purpose. In [START_REF] Meier | Analysis of pseudo random sequences generated by cellular automata[END_REF], Meier and Staffelbach proposed a correlation attack to reverse the PRS generated by rule 30 although it passes classical statistical tests like the ones proposed in [START_REF] Knuth | Seminumerical Algorithms[END_REF].

More recently, in [START_REF] Martin | A Walsh exploration of elementary CA rules[END_REF], we have used a Walsh transform to explore the set of the 256 elementary rules. The Walsh transform is a well-known tool in the field of cryptology for studying the correlation-immunity of Boolean functions: Xiao and Massey [START_REF] Xiao | A spectral characterization of correlation-immune combining functions[END_REF] have characterized the notion of correlation-immunity with the Walsh transform. We have applied this technique to the pseudo-random sequences generated by all of the 256 binary rules and we provide evidence that there does not exist a non-linear rule which generates a correlation-immune pseudo-random sequence. Thus, we state Theorem 6. Theorem 6. [START_REF] Martin | A Walsh exploration of elementary CA rules[END_REF] There is no non-linear correlation-immune elementary CA.

And, according to Theorem 2, we can state that: Corollary 3. There is no elementary CA which can serve as PRS generator.

So, does Theorem 6 annihilate any hope to design a good PRG by the means of CA? Not necessarily. Next section recalls the approach initiated by Tomassini and Sipper and section 5 describe another way of generating PRS with LHCA.

PRS generation by HCA

The cellular programming approach

Tomassini and Sipper [START_REF] Sipper | Co-evolving parallel random number generators[END_REF] proposed to use HCA for generating better PRS. In this model, the rules are obtained by an evolutionary approach (a genetic algorithm). They have designed a cellular programming algorithm for cellular automata to perform computations, and have applied it to the evolution of pseudo-random sequence generators. Their genetic algorithm uses Koza's entropy E h = -k h j=1 p hj log 2 p hj where k denotes the number of possible values per sequence position, h a subsequence length and p hj is a measured probability of occurrence of a sequence h j in a pseudo-random sequence. It measures the entropy for the set of k h probabilities of the k h possible subsequences of length h. The entropy achieves its maximal value E h = h when the probabilities of the k h possible sequences of length h are all equal to 1/ℓ h , where ℓ h denotes a number of possible states of each sequence. They have selected four rules of radius 1 for use in non-uniform cellular automata. The best rules selected by the genetic algorithm were rules 90, 105, 150 and 165 (which are all linear, a clear drawback).

A series of tests (including χ 2 test, serial correlation coefficient, entropy and Monte Carlo, but no correlation-immunity analysis) were made with good results, showing that co-evolving generators are at least as good as the best available CA randomizer. The authors also use elementary rules which we proved to be not correlation-immune. This was further investigated in [START_REF] Seredynski | Cellular automata computations and secret key cryptography[END_REF].

Following the same kind of approach, Seredynski et al. in [START_REF] Seredynski | Cellular automata computations and secret key cryptography[END_REF] have generalized the selection process to radius 2 rules. They use then both radius 1 and radius 2 rules in hybrid cellular automata. The rules selected by their genetic algorithm were 30, 86, 101 and 869020563, 1047380370, 1436194405, 1436965290, 1705400746, 1815843780, 2084275140 and 2592765285.

Their new set of rules was tested by a number of statistical tests required by the FIPS 140-2 standard [START_REF]Security requirements for cryptographic modules[END_REF] but no correlation-immunity analysis was made either.

The synthesis approach

This approach follows the synthesis algorithm proposed in [START_REF] Cattell | Synthesis of one-dimensional linear hybrid cellular automata[END_REF]. They propose a method for the synthesis of a HCA from a given irreducible polynomial over F 2 . The same problem for LFSR is well known as it can be directly obtained from the transition matrix. Furthermore, there is a one to one correspondence between LFSR's and polynomials. For CA, in general, a characteristic polynomial is not sufficient to uniquely determine the CA from which it was computed.

If we consider the characteristic polynomial ∆ of the HCA (assumed to be irreducible), with α a root in F 2 n . All n roots of ∆ lie in F 2 n . The roots α, α 2 , α 2 2 , . . . , α 2 n-1 are distinct and ∆ can be factored in

F 2 n as (x -α)(x -α 2 )(x - α 2 2 ) . . . (x -α 2 n-1 ).
Product of irreducibles Given p and q two irreducible characteristic polynomials and P and Q their respective transition matrix, on can build the transition matrix corresponding to p • q. It can be defined by blocks as: P 0 0 Q . This operation corresponds to the concatenation of LHCA [START_REF] Cattell | One-dimensional LHCA: their synthesis, properties and applications in VLSI testing[END_REF]. They quoted that it permits to concatenate primitive machines for forming machines of much longer lengths.

Application: boolean functions evaluation

There is a well-known dictionary between, on the one hand, boolean functions in n variables, and binary sequences of period 2 n . More specifically, if f is such a function, and if we denote by i the base 2 expansion of i we can define a sequence by the rule

s f (i) = f (i), for i ≤ 2 n -1.
Many interesting boolean functions can be cast under the form f (x) = T r(ax + bx s ), where a, b are scalars of the extension field F 2 n and T r the trace function from F 2 n down to F 2 . In the case where n is odd and the Walsh Hadamard transform takes only three values they are the so-called plateaued boolean function of order n -1 [START_REF] Zheng | On plateaued functions[END_REF] also known as almost optimal or semi-bent.

They are the traces of so-called almost bent AB functions [START_REF] Carlet | Codes, bent functions and permutations suitable for DES-like cryptosystems[END_REF]. For monomials AB functions, the most famous exponents s are in the conjecturally exhaustive list of Gold, Kasami, Welch, Niho (see Table 1). In all these cases, an upshot of the theory of Name s Condition Gold

2 i + 1 i ∧ m = 1, 1 ≤ i ≤ m/2 Kasami 2 2i -2 i + 1 i ∧ m = 1, 1 ≤ i ≤ m/2
Welch 2 (m-1)/2 + 3 Niho 2 2r + 2 r -1 r = t/2 for t even r = (3t + 1)/2 for t odd with 1 ≤ r ≤ m = 2t + 1 Table 1. Exponents of AB monomials. Mattson-Solomon polynomials [10, p.249] is that the parity check polynomials of the attached cyclic codes (or, essentially, the connection polynomial of the LFSR) is of the form m α m α s , where α generates F 2 n over F 2 . A fast algorithm to compute the minimal polynomials of elements in finite field extensions is given in [START_REF] Augot | On the computation of minimal polynomials, cyclic vectors, and frobenius forms[END_REF].

Conclusion

We have used the synthesis approach to give an effective CA-realization of classical pseudo-random sequences of cryptographic quality. The main interest of this work would be to give an hardware implementation. The target hardware model of CAs is the Field Programmable Gate Arrays (known as FPGAs). FPGAs are now a popular implementation style for digital logic systems and subsystems. These devices consist of an array of uncommitted logic gates whose function and interconnection is determined by downloading information to the device. When the programming configuration is held in static RAM, the logic function implemented by those FPGAs can be dynamically reconfigured in fractions of a second by rewriting the configuration memory contents. Thus, the use of FPGAs can speed up the computation done by the cellular automata. Putting all together allows high-rate pseudo-random generation of good quality.

Fig. 2 .

 2 Fig. 2. Evolution of CA30 on a ring with N = 8 cells.

⋆ This work was supported by the french ANR program NUGET.


Fabio Guttler
email: fabio.guttler@teledetection.fr

Dino Ienco
email: dino.ienco@teledetection.fr

Maguelonne Teisseire
email: maguelonne.teisseire@teledetection.fr

Jordi Nin
email: nin@ac.upc.edu

Pascal Poncelet
email: pascal.poncelet@lirmm.fr

Towards the Use of Sequential Patterns for Detection and Characterization of Natural and Agricultural Areas

Keywords:

Nowadays, a huge amount of high resolution satellite images are freely available. Such images allow researchers in environmental sciences to study the di↵erent natural habitats and farming practices in a remote way. However, satellite images content strongly depends on the season of the acquisition. Due to the periodicity of natural and agricultural dynamics throughout seasons, sequential patterns arise as a new opportunity to model the behaviour of these environments. In this paper, we describe some preliminary results obtained with a new framework for studying spatiotemporal evolutions over natural and agricultural areas using k-partite graphs and sequential patterns extracted from segmented Landsat images.

Introduction

Several regions over the earth are composed of complex landscapes with regard to land cover and land use. Outside urban areas, we commonly observe mosaics of natural, semi-natural and agricultural areas. In Europe, mapping and monitoring those areas became a major issue and several procedures have been established for guiding and controlling such tasks. With regard to natural areas we can cite the Habitats Directive (92/43/EEC) and the associated Natura 2000 network of protected sites. In this context, member states must report the conservation status of the habitats within their territory every six years [START_REF] Vanden Borre | Integrating remote sensing in natura 2000 habitat monitoring: Prospects on the way forward[END_REF]. Concerning agriculture, Land Parcel Identification Systems (LPIS) emerged since 1992 (Council Reg. No 3508/1992). LPIS is used as a reference for annual declaration by farmers, for administrative and cross-checks, and on the spot controls (including control with remote sensing) [START_REF] Trojáček | Detailed mapping of agricultural plots using satellite images and aerial orthophoto maps[END_REF].

The detection of natural and agricultural areas from remote sensing data is a key point for land cover mapping. Image classification is a widespread method for mapping land cover. The overall objective is to categorize all pixels in an image into land cover classes [START_REF] Lillesand | Remote Sens. and Image Interpret[END_REF]. Usually, image classification is based on the di↵erential spectral response of land surfaces (radiance values recorded at pixel level). However, most of classifications are based on single image and present lot of bias in their results, or require extensive ground truth data in order to attempt a high accuracy.

Temporal pattern recognition can furnish complementary information to feature identification. Actually, natural vegetation and agricultural crop present distinct evolutions during phenological cycles and growing seasons. In that light, multi-date imagery can enhance the pertinence and the accuracy of land cover detection. Time series analysis comprises methods for analysing temporal data [START_REF] Petitjean | Spatio-temporal reasoning for the classification of satellite image time series[END_REF], several images in our case, in order to extract meaningful statistics and other characteristics of the data. Temporal data allows researchers to create a data model for analysing past values and forecasting future ones [START_REF] Nin | Towards the evaluation of time series protection methods[END_REF]. Specifically, this work is focus on natural and agricultural areas modelling over time. Sequential pattern mining is the part of time series analysis concerned with finding statistically relevant patterns between data examples where the values are delivered in a di↵erent time moments [START_REF] Nin | Speed up gradual rule mining from stream data! a b-tree and owa-based approach[END_REF].

Optical remote sensing, such as Landsat images are commonly used in environmental researches. Image processing techniques are usually grouped in Pixel based analysis or Object based image analysis (OBIA). Several time series pixelbased approaches have been proposed, but OBIA studies rarely uses multitemporal data [START_REF] Petitjean | Spatio-temporal reasoning for the classification of satellite image time series[END_REF].

In this paper we combine OBIA with sequential pattern mining to create a k-partite graph to represent natural and agricultural areas. Our approach starts segmenting multi-date satellite images over a same area. From these segmentations, a k-partite graph is built considering objects as nodes and representing object image overlap as edges. Later, in order to determine the optimum bounding box of each object and study its dynamics, we propose the creation of k-partite subgraphs. Finally, subgraphs are supplied to an expert in order to be categorized.

Paper Organisation

This paper is organised as follows: Firstly, we introduce some concepts about optical remote sensing in Section 2. Section 3 describes k-partite graphs and sequential pattern representation. Our technique to maximise image coverage is introduced in Section 4. Experimental findings about agricultural and natural areas are depicted in Section 5. Section 6 concludes and draws future works.

Remote Sensing Satellite Images

In general we talk about remote sensing when the acquisition of data is done without making any physical contact. This is the case of Earth observation plat-forms such as Landsat. The optical sensors onboard Landsat acquires multispectral images (composed by multiple bands), each band representing a portion of the electro-magnetic spectrum. For this study, we used Landsat Thematic Mapper (TM) images without the thermal band. Remote sensing images are usually characterized by di↵erent resolutions and other technical characteristics, in our case they are as follows:

-Spatial resolution. The size of a pixel in a raster image is 30x30 meters.

-Spectral resolution. The spectral bands are six: blue, green, red, near infrared (NIR) and two parts of short wavelength infrared (SWIR-1 and 2). -Radiometric resolution. The sensor is able to distinguish is 256 intensities of radiation (8 bits). -Swath width. The scene size is a ground square of about 185 x 185 km.

Apart from these parameters another important issue in time-series studies is the frequency of revisits by the satellite. The revisit time of Landsat is of 16 days, but in practice the number of useful images is lower, mostly due to unsuitable whether conditions during satellite acquisition. For this study, we selected six images covering of the same geographic area (the Natura 2000 Lower Aude Valley site, located in south of France) between February and September 2009.

To avoid errors and misplacement over time, all images were already georeferenced, as well as radiometrically and atmospherically corrected by CESBIO-CNES [START_REF] Hagolle | Séries temporelles de produits landsat de niveau 2a: Manuel de l'utilisateur des données produites au cesbio[END_REF]. Additional fine spatial positioning corrections were applied in order to keep the spatial shift between all the time stamps less than one pixel.

Segmentation

Image segmentation is a fundamental step in OBIA and its consists in merging pixels into objects clusters [START_REF] Baatz | Progressing from object-based to objectoriented image analysis[END_REF]. Objects (or segments) are regions generated by one or more criteria of homogeneity in one or more dimensions of a feature space [START_REF] Blaschke | Object based image analysis for remote sensing[END_REF]. The aim of segmentation is to create a new representation of the image more meaningful and easier to analyse. This approach is similar to human visual interpretation of digital images, which works at multiple scales and uses colour, shape, size, texture, pattern and context information [START_REF] Lillesand | Remote Sens. and Image Interpret[END_REF].

Image segmentation results in a set of objects that collectively cover the entire image without any overlapping. With respect to the homogeneity criteria, adjacent objects are expected to be significantly di↵erent between them. In this work we use the software eCognition Developer 8.8.1 for image segmentation (multiresolution segmentation algorithm). Only the pixels within the boundaries of the Lower Aude Valley Natura 2000 site (4,842 ha) were used for segmentation and further processing steps. Nine raster layers (radiometric or 'colour' information) were used simultaneously for image segmentation. Six of them correspond to the Landsat spectral bands and the other are spectral indices. Spectral indices are commonly used in remote sensing as they can be helpful for detecting and characterizing some specific features, like vegetation, soil, water, etc. We used the Normalized Di↵erence Vegetation Index (NDVI) [START_REF] Rouse | Monitoring vegetation systems in the great plains with erts[END_REF], the Normalized Di↵erence Water Index (NDWI) [START_REF] Jackson | Vegetation water content mapping using landsat data derived normalized di↵erence water index for corn and soybeans[END_REF] and the Visible and Shortwave Infrared Drought Index (VSDI) [START_REF] Zhang | Vsdi: a visible and shortwave infrared drought index for monitoring soil and vegetation moisture based on optical remote sensing[END_REF]. In order to obtain objects of interest related to natural and agricultural areas, we conceived a segmentation rule-set composed of 3 main steps:

1. Medium-coarse segmentation to delineate general zones (colour and shape components combined but "color > shape"): about 170-200 objects 2. Very fine segmentation focused on colour component: about 6,000 objects 3. Medium-fine segmentation with balanced weights for colour and shape components: about 500-600 objects This process was independently applied for each Landsat image. The last segmentation level (3) was then exported from each time-stamp and used as input for the subsequent processing steps. Figure 1 illustrates the segmentation layer obtained for the time stamp of 10 July 2009.

Sequential Patterns and k-partite Graphs

This section describes the traditional sequential pattern mining problem and high-lights the need for a specific way to handle remote sensing temporal information by using k-partite graphs.

k-partite Graphs

A k-partite graph is a graph G = (V, E) with vertex set V and edge set E, whose graph vertices V can be partitioned into k disjoint sets V k so that no two vertices within the same set are adjacent. In this work, we assume that any edge e ij has a weigh equal to w ij where i and j correspond to two vertexes v i 2 V i and v j 2 V j in two consecutive sets (layers).

Sequential Patterns

Sequential patterns were introduced in [START_REF] Agrawal | Mining sequential patterns[END_REF] and in general, they are considered as an extension of the concept of frequent itemset [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF] having timestamps associated to items. Sequential pattern mining aims at extracting sets of items commonly associated over time.

The problem of generating and mining sequential patterns as a n-bipartite graph is defined as follows.

Let X be a set of distinct items. An itemset is a subset of items, denoted by I = (i 1 , i 2 , . . . , i n ), for 1  j  n, i j 2 X. A k-partite graph is built as an ordered list of itemsets, denoted by hI 1 , I 2 , . . . , I k i, where I i 2 X for 1  i  n and correspond with each V i nodes set.

Once graph nodes are defined, it is required to define the edges and their corresponding weights. To do that, a function f (I i , I j ) = w ij is defined, if f (I i , I j ) < 0 it means that there is one edge between a pair of nodes v i and V j with weight equal to w ij .

Representation of natural and agricultural areas

In order to represent natural and agricultural area evolutions over time, it is proposed to create a sequential pattern by using a k-partite graph representation in the following way:

1. Segmented images are the di↵erent itemsets hI 1 , I 2 , . . . , I k i, note that, the image timestamp provide us the temporal information for creating the different nodes sets and therefore the k-partite graph. 2. Objects extracted from the segmented images allow us to generate the items i l of each itemset I j . 3. To create the edges, the function f

(I i , I j ) = |p i l \ p i l+1 |/|p i l [ p i l+1 | is cal- culated,
where p i l stands for the list of pixels of the item i l . Within this global k-partite graph, it is possible to extract a concrete area of interest evolution creating a sub-graph chosen one or more items. All this process is illustrated in Figure 2.

Coverage Algorithm

In this paper we have implemented the following algorithm to (completely) cover one region of interest with sequential patterns: firstly, we detect the previously unprocessed areas of the image where frequent changes in the pixel values are produced. After that, we select from the first image the objects included in the previously selected image area. Later, we generate a sequential pattern and mark this part of the image as processed and recompute object variability. The algorithm stops when the region of interest has been completely processed.

Specifically, we implemented the above coverage algorithm as it is described in Algorithm 1. Firstly, we create a bounding box list containing all the objects of all the images (lines 2-5), at the same time, we compute the variability of each bounding box (object). Bounding box importance is the amount of other objects overlapping with at least one pixel with the current bounding box, this computation is described in Algorithm 2.

Once the bounding box list is created, all objects are ranked with regards their importance. The importance of a bounding box corresponds with the number of objects in any layer overlapping with at least one pixel with the bounding box. After that, the most important bounding box is selected (line 7) and recover the objects overlapping with such bounding box in the first image (line 8). Then, a new pattern is generated as it is described in Section 3.2 (line 9). After that, we update bounding boxes importance, subtracting the processed objects from the importance of each bounding box, if bounding box importance is equal to 0, it is removed from the list (Algorithm 3). This loop is repeated until the bounding box list is empty.

T 0 T 1 T 2 T 3 T 4 T 5

Experiments

In order to evaluate the pertinence of our approach, graphical results were analyzed by a remote sensing expert with field knowledge of the study area. For this analysis, all subgraphs (331 in total) were plotted with an associated unique ID allowing to locate each bounding box in an Geographic Information System (GIS). The expert preselected (without taking into account subgraph results) about 20 regions of interest (ROI) containing both natural and agricultural areas. Then, the subgraphs corresponding to this set of ROI were identified and analyzed. The following subsection details the main di↵erences between the graphs representing spatio-temporal evolutions of agricultural and natural areas within our study area. Agricultural area. In general, agricultural areas present less complex subgraphs as the fields are usually rectangular and the boundaries of the objects remain very similar from one time stamp to the next one. Even if subgraph structure is quite simple (i.e. Figure 4 -subgraph 1), radiometric information varies sharply throughout the time series. Subgraph 1 (Figure 4) illustrates the temporal pattern of a winter culture (cereal crop in this example). The parcel presents an important augmentation of the NDVI from late winter up to spring, which corresponds to the growing season. The crop is harvest in late spring or early summer, generating a brutal change in the spectral response (between T2 and T3). Afterwards, the field is characterized by bare soil, which remains globally stable until the end of the time series.

Pattern analysis examples

Natural area. In natural areas the boundaries are usually much smoother than in agricultural fields and the objects rarely have rectangular shapes. As a consequence, graph structure tends to be more complex and can englobe several objects per layer. Such kind of graph can also show some drastic radiometric changes, but in general it concerns only a small portion of the area covered by the bounding box. In other words, natural habitats are characterized by more gradual and smoothed evolutions. Subgraph 2 of figure 4 illustrates the temporal pattern of a temporary lake. Usually flooded during winter and spring, this lake becomes progressively dry starting from early summer. Aquatic natural habitat is gradually replaced by some pioneer terrestrial communities, dominated by annual Salicornia sp. or Suaeda sp. Depending mainly on the local variations of soil salinity, the vegetation will cover the former flooded area in di↵erent ways.

First developments of such particular habitats can be detected since T2 and they should subsist up to the next flood event (usually during late autumn or winter).

Discussion

Throughout one year, spatio-temporal evolutions of agricultural and natural areas are enough dissimilar and it is possible to an expert to separate them by analysing the time series of satellite images. In a more automatic way, the subgraphs generated by our approach synthesizes complex spatio-temporal evolutions and can be useful for this kind of time-consuming task. Actually, subgraph structure associated to temporal behaviour of object's radiometry provides complementary and pertinent information allowing detailed analysis. At this stage, our coverage algorithm is not able to cover 100% of the study area throughout the time series. Most of times, small objects are not included in the subgraphs. However, this issue should be improved as the spatial gaps concerns sometimes also medium-size objects representing interesting natural or agricultural areas. Another point to improve is spatial redundancy. We verified that some subgraphs presents high rates of similar objects, what means that the same area is covered by more than one subgraph. Improving the bounding box selection should reduce such spatial redundancy.

Conclusions

In this paper we have described a complete framework for studying evolving natural and agricultural areas using satellite images information and k-partite graph sequential patterns. We have shown that using our approach is possible to cover a great part of the study area (93.56%) and to analyze in detail a concrete region of interest. We have also verified with one domain expert that the obtained sequential patterns are meaningful. As a future work, we would like to develop some clustering algorithms for k-partite graphs to help the expert to post-process the results providing some groups of patterns with similar behaviour, instead of individuals patterns without relation among them.

Fig. 1 .

 1 Fig. 1. Segmentation example representing the time stamp of 10 July 2009.

Fig. 2 .r e a c h layer l 2 g do 4 f o r e a c h node n 2 l do 5

 25 Fig. 2. Overall Sequential Patterns Extraction Process.

Algorithm 2 :r e a c h layer l 2 g do 4 f o r e a c h node n 2 r e a c h boundingBox b 2 bb do 3 f o r e a c h layer l 2 p do 4 f o r e a c h node n 2

 222 Importance Data: g: main graph, b:bounding box Result: l do /* if one pixel of node n l do /* if one pixel of node n

Figure 3

 3 Figure 3 shows the obtained coverage within the perimeter of our study area. Each polygon represents the spatial extent of a specific bounding box, what means that the associated pattern has a bigger spatial coverage taking into account all temporal layers. Polygon color indicates from which time stamp the bounding box was extracted (as assigned on the legend). In total 331 bounding box were selected (about 50 -60 per timestamp), it corresponds to an overall coverage of 93.56% of the study area.

Fig. 3 .

 3 Fig. 3. Obtained spatial coverage.

Fig. 4 .

 4 Fig. 4. Temporal Pattern examples: Agricultural area (subgraph 1) and Natural area (subgraph 2)

Acknowledgments. This work was supported by public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the program "Investissements d'Avenir" managed by the French National Research Agency. This work was also partially supported by the Ministry of Science and Technology of Spain under contract TIN2012-34557 and by the BSC-CNS Severo Ochoa program (SEV-2011-00067).