L Albera

S Kitić

N Bertin

G Puy

Rémi Gribonval

BRAIN SOURCE LOCALIZATION USING A PHYSICS-DRIVEN STRUCTURED COSPARSE REPRESENTATION OF EEG SIGNALS

Keywords: Brain source localization, EEG, cosparsity, synchronous current activities

Localizing several potentially synchronous brain activities with low signal-to-noise ratio from ElectroEncephaloGraphic (EEG) recordings is a challenging problem. In this paper we propose a novel source localization method, named CoRE, which uses a Cosparse Representation of EEG signals. The underlying analysis operator is derived from physical laws satisfied by EEG signals, and more particularly from Poisson's equation. In addition, we show how physiological constraints on sources, leading to a given space support and fixed orientations for current dipoles, can be taken into account in the optimization scheme. Computer results, aiming at showing the feasability of the CoRE technique, illustrate its superiority in terms of estimation accuracy over dictionarybased sparse methods and subspace approaches.

INTRODUCTION

Electrical potentials produced by neuronal activity can be measured at the surface of the head using ElectroEncephaloGraphy (EEG). The sources of this neuronal activity (during either cognitive or pathological processes) can be localized provided that a model of sources and a model of volume conductor are defined. Most source localization algorithms use one of the two following source models: the point source model, which explains the data with a small number of equivalent current dipoles and the distributed source model, which uses thousands of dipoles. Whereas the latter allows for an estimation of the spatial extent of the source, it requires to make assumptions about the spatial source distribution, which may lead to blurred (or even distorded) solutions [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF]. On the other hand, the former often gives helpful first approximations and superior performance in environments where there are few sources which are clustered [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF]. Regarding head models, The work of S. Kitić, G. Puy and R. Gribonval was funded by the FP7 European Research Council Programme, PLEASE project, under grant ERC-StG-2011-277906. they aim at representing geometrical and electrical properties of the different tissues composing the volume conductor. Various models were proposed going from concentric homogeneous spheres with isotropic conductivities to realistically shaped models with refined tissue conductivity values [START_REF] Vatta | Realistic and spherical head modeling for EEG forward problem solution: A comparative cortex-based analysis[END_REF].

Numerous methods were developed to localize equivalent current dipoles from EEG recordings. Among them, beamforming techniques [START_REF] Piotrowski | Reduced-rank neural activity index for eeg/meg multi-source localization[END_REF], subspace approaches such as MU-SIC (MUltiple-SIgnal Classification) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] and sparse methods [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF] are the most popular. One can mention the RapMUSIC (Recursively applied MUSIC) [START_REF] Mosher | Source localization using Recursively Applied and Projected (RAP) music[END_REF] and FO-D-MUSIC (Fourth Order Deflationary MUSIC) [START_REF] Albera | Brain source localization using a fourth-order deflation scheme[END_REF] algorithms, which are sequential versions of MUSIC based on Second Order (SO) and Fourth Order (FO) statistics, respectively. Regarding the dictionary-based sparse techniques, the most famous is MCE (Minimum Current Estimate) [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF], which computes minimum ℓ 1 -norm estimates. Note that cosparse approaches were recently proposed applying for instance a discrete Gabor transform [START_REF] Gramfort | Timefrequency mixed-norm estimates: Sparse m/eeg imaging with non-stationary source activations[END_REF] as analysis operator [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] to the current source activities.

In this paper, we show how to derive an analysis operator derived from physical laws satisfied by EEG signals, and more particularly from Poisson's equation. As a result, we propose a new dipole localization method, named CoRE (physics-driven structured Cosparse Representation of EEG signals), taking into account physiological constraints on sources leading to a given space support and fixed orientations. Computer simulations demonstrate the feasability and the performance of the CoRE approach in comparison to the MCE technique especially for low Signal-to-Noise Ratio (SNR) values. We also illustrate the ability of CoRE in dealing with synchronous brain sources and overcoming the RapMUSIC and FO-D-MUSIC algorithms in such a practical context.

ASSUMPTIONS AND PROBLEM FORMULATION

It is commonly admitted that the electrical potential v(r) recorded at location r of the head mostly reflects the activity of pyramidal cells located in the gray matter and oriented perpendicularly to the cortical surface. This activity is gen-erally modeled by current dipoles. Given the geometry and the scalar field {σ(r)} of electrical conductivities within the head, Poisson's equation [START_REF] Mohr | Comparing iterative solvers for linear systems associated with the finite difference discretisation of the forward problem in electroencephalographic source analysis[END_REF] allows us to establish a link between the electrical potential v(r) and the electrical activity of the current dipoles:

div r (σ(r) ∇ r v(r)) = ρ∈S s(ρ)(δ(r -ρ -) -δ(r -ρ + ))/d (1)
where ρ -and ρ + denote the positions of the two monopoles representing the current sink and source, respectively, d stands for the distance between the two monopoles, s(ρ) is the electrical activity of the dipole located at position ρ and S denotes the gray matter. The location parameter ρ of the considered dipole is typically chosen half way between the two monopoles. Let's consider the following assumptions: Then, the brain source localization problem consists in identifying the locations ρ p of the P dominant current dipoles belonging to S from the N -dimensional vector x = [v(r 1 ), . . . , v(r N )] T .

A1. N potential values v(r n ) recorded

PHYSICS-DRIVEN STRUCTURED COSPARSE REPRESENTATION OF EEG SIGNALS

Now let's see how the brain source localization problem can be reformulated as a cosparse analysis model fitting problem based on a physics-driven structured cosparse representation of EEG signals.

Discretization of Poisson's equation

The CoRE approach described in section 4 requires an appropriate discretization of Poisson's equation in order to derive an equation of the form Ω v = z from [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF] where Ω is a matrix denoting the so-called linear analysis operator [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] and where z is a sparse vector. The Finite Difference Method (FDM) and Finite Element Method (FEM) can be used to perform such a discretization. Note that these two methods allow us to take realistic head models into account such as those obtained from anatomical imaging modalities: Computed Tomography (CT) and structural Magnetic Resonance Imaging (sMRI) for instance.

As an example, and without loss of generality, we consider the FDM in the following, and more particularly the FDM implemented by Witwer et al. [START_REF] Witwer | The effect of media inhomogeneities upon intracranial electrical fields[END_REF] applying Kirchhoffs law at each node of the cubic grid. It is noteworthy that Witwer et al. gave an easy way of taking into account the skull-air boundary. The nodal points (voxels) exterior to the insulated boundary are assigned a conductivity that is so small that the flux leakage across the boundary can be made as small as desired. Another main advantage of Witwer's FDM is the way it is treating the internal boundaries and inhomogeneities, say by simply varying the conductivity values. Discretizing the volume as a (L × L × L) uniform three-dimensional grid with spacing h, and applying Witwer's FDM to equation ( 1), we can express the total current flow value at a given voxel as a linear combination of the potentials measured in a close neighbourhood of this voxel (see equation ( 17) of appendix).

Toward a cosparse analysis model

The link between the total current flow value and the potentials derived from the discretization of Poisson's equation can be easily written in the matrix form A v = i where v and i are two vectors containing the potential and total current flow values, respectively, at the different nodes of the threedimensional grid. The (L 3 × L 3 ) matrix A is symmetric, positive semidefinite, rank deficient by one and sparse with only seven non-zero components in each row [START_REF] Mohr | Comparing iterative solvers for linear systems associated with the finite difference discretisation of the forward problem in electroencephalographic source analysis[END_REF]. Generally, instead of considering the singular linear system A v = i, another possibility is to transform it into a regular one and solve this instead. The regular system is chosen such that its unique solution belongs to the set of solutions of the original singular system. As described in [START_REF] Mohr | Comparing iterative solvers for linear systems associated with the finite difference discretisation of the forward problem in electroencephalographic source analysis[END_REF], the easiest approach is to fix the value of the potential to zero in one voxel. The special structure of the matrix A then allows us to cancel the corresponding row and column in A and also the respective entry in the right-hand side vector i. This leads to a system for which the (L 3 -1 × L 3 -1) resulting matrix is symmetric and positive definite. By slight abuse of notation we continue to name this matrix A. The solution of this system obviously solves the initial system with a zero potential value in the respective voxel.

Let ℓ -and ℓ + index the elements of i giving the amplitude of the monopoles located at positions ρ -and ρ + , respectively, of the grid. According to assumptions (A2), (A3) and (A4), and knowing the geometry of the cortex, it is possible to build a set of Q different couples (ℓ - q , ℓ + q ) indexing the Q dipoles of S. Indeed, by covering the surface of the cortex with Q monopoles indexed by the integers ℓ - q , we can deduce the Q corresponding integers ℓ + q using (A4). We can thus factorize vector i as i = B z where B = (B ℓ1,ℓ2 ) is a (L 3 -1 × L 3 -1) sparse matrix defined by:

B ℓ1,ℓ2 =    1 if ℓ 1 = ℓ 2 -1 if ℓ 1 = ℓ + q and ℓ 2 = ℓ - q 0 otherwise (2)
and where z is a (L 3 -1)-dimensional P -sparse vector with L 3 -1-Q known zero elements: the ℓ-th element of z is zero if ℓ = ℓ - q for q ∈ {1, . . . , Q}. Consequently, matrix B conveys our knowledge about the orientation of the Q dipoles of S and the amplitude of the monopoles associated with these dipoles. It is also an invertible matrix, since it is possible to permute its rows and columns such that a lower triangular matrix with ones on the diagonal is obtained. Regarding vector z, its non-zero elements represent the amplitude of monopoles restricted to the cortical surface.

The brain source localization defined in section 2 can then be reformulated as a cosparse analysis model fitting problem given by:

Ω v = z x = M v (3)
where the analysis operator Ω is given by Ω = B -1 A and the matrix M is a (N × L -1) row-reduced identity matrix.

THE CoRE ALGORITHM

The purpose of this section is to show how the cosparse analysis model fitting problem presented in section 3 can be solved, giving birth to the CoRE algorithm. First, we reformulate the cosparse analysis model fitting problem (3) as the following convex problem:

minimize v ||Ω 1 v|| 1 + λ||Ω 2 v|| 2 2 (4) subject to M v = x.
where Ω 1 is the (Q × L -1) submatrix of Ω obtained by extracting the rows of Ω corresponding to the support set S, whereas Ω 2 corresponds to the rows indicated by the complementary set S. By choosing the appropriate weight λ, the cosparse solution of the optimization problem (4) will fulfill the assumptions (A1) to (A5). Namely, ||Ω 1 v|| 1 will promote sparsity at the surface of the cortex, while λ||Ω 2 v|| 2 2 will attenuate the signal in the other regions. The linear constraints M v = x assure that the model fits the electrode measurements. Depending on the resolution of the cubic grid tuned by L, the problem can reach considerably large scale. Therefore, we propose an efficient first order method in order to solve it, through the Alternating Directions Method of Multiplier (ADMM) framework [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. First, we reformulate the problem (4) as follows:

minimize v,z1,z2 z 1 1 + λ z 2 2 2 + f (v) (5)
subject to

Ω 1 Ω 2 v = z 1 z 2
where f denotes an indicator function of the set

C = {v ∈ | L-1 | M v = x} given by f (v) = 0 for v ∈ C and f (v) = +∞ otherwise.
Then, the augmented Lagrangian has the following form:

L ρ1,ρ2 (v, z 1 , z 2 , µ 1 , µ 2 ) = z 1 1 + λ z 2 2 2 + f (v) + µ 1 T (Ω 1 v -z 1 ) + ρ 1 2 Ω 1 v -z 1 2 2 + µ 2 T (Ω 2 v -z 2 ) + ρ 2 2 Ω 2 v -z 2 2 2 , (6)
or after substitution u 1 = 1 ρ1 µ 1 and u 2 = 1 ρ2 µ 2 , the following convenient scaled form:

L ρ1,ρ2 (v, z 1 , z 2 , u 1 , u 2 ) = z 1 1 + λ z 2 2 2 + f (v) + ρ 1 2 Ω 1 v -z 1 + u 1 2 2 - ρ 1 2 u 1 2 2 + ρ 2 2 Ω 2 v -z 2 + u 2 2 2 - ρ 2 2 u 2 2 2 . ( 7
)
We now perform an iterative alternating minimization with respect to each of the primal variables v, z 1 and z 2 . The minimization with respect to v requires to solve an equality constrained least squares problem:

v (k+1) = arg min v L ρ1,ρ2 (v, z (k) 1 , z (k) 2 , u (k) 1 , u (k) 2 ) = arg min v Hv -w 2 2 subject to M v = x. (8)
where

H = [ √ ρ 1 Ω 1 T , √ ρ 2 Ω 2 T ] T and w = [ √ ρ 1 (z (k) 1 - u (k) 1 ) T , √ ρ 2 (z (k) 2 -u (k) 2 ) T ] T . Here v (k) , z (k) 1 , z (k) 2 , u (k) 1 and u (k) 2
stand for the estimates of the solutions minimizing [START_REF] Albera | Brain source localization using a fourth-order deflation scheme[END_REF], computed at iteration k. Solving the problem ( 8) is straightforward when the solution is expressed as the following sum [13, section 16.2]:

v (k+1) = v p + Kv Null (9)
where v p is a particular solution of the underdetermined constraint system M v = x and v Null is the null space component with K a null space basis of M . Due to the particular structure of M , which is a row-reduced identity matrix, the particular solution is simplified to v p = M T x and the columns of K can be simply computed by transposing the rows of the (L -1 × L -1) identity matrix, which do not appear in M . Then, by variable substitution in the objective, we obtain:

v Null = arg min y H K y + H M T x -w 2 2 . ( 10
)
Updating the z

(k)
1 variable is equivalent to computing the proximity operator for the ℓ 1 -norm [START_REF] Bach | Convex optimization with sparsityinducing norms[END_REF]:

z (k+1) 1 = arg min ż L ρ1,ρ2 (v (k+1) , z 1 , z (k) 2 , u (k) 1 , u (k) 2 ) = arg min ż z 1 1 + ρ 1 2 z 1 -Ω 1 v (k+1) -u (k) 1 2 2 z (k+1) 1 = prox 1 ρ 1 ℓ1 Ω 1 v (k+1) + u (k) 1 , (11)
which is an element-wise soft thresholding (see [12, section 4.4.3] for more details). The third update concerns the minimization of a weighted sum of two quadratic functionals:

z (k+1) 2 = arg min z2 L ρ1,ρ2 (v (k+1) , z (k+1) 1 , z 2 , u (k) 1 , u (k) 2 ) = arg min z2 λ z 2 2 2 + ρ 2 2 z 2 -Ω 2 v (k+1) -u (k) 2 2 2 z (k+1) 2 = ρ 2 (2λ + ρ 2 ) Ω 2 v (k+1) + u (k) 2 . ( 12
)
The dual variables are updated in the final stage using the following rules:

u (k+1) 1 = u (k) 1 + Ω 1 v (k+1) -z (k+1) 1 u (k+1) 2 = u (k) 2 + Ω 2 v (k+1) -z (k+1) 2 (13)
This is the standard way to compute the updates for the scaled Lagrangian multipliers (check [12, section 3.1.1]).

The stopping criterion is based on the primal and dual residuals:

r (k+1) pri = Ω 1 Ω 2 v (k+1) - z (k+1) 1 -u (k+1) 1 z (k+1) 2 -u (k+1) 2 r (k+1) dual = ρ 1 Ω 1 T z (k+1) 1 -z (k) 1 ρ 2 Ω 2 T z (k+1) 2 -z (k) 2 , (14)
with appropriate thresholds:

ǫ (k+1) pri = Qǫ A + ǫ R max{ Ω v (k+1) , z (k+1) 1 + z (k+1) 2 } ǫ (k+1) dual
= √ L -1ǫ A + ǫ R ( ρ 1 Ω 1 T u (k+1) 1 2 + ρ 2 Ω 2 T u (k+1) 2
2 ).

where ǫ A and ǫ R denote the absolute and relative tolerances, which we empirically set to 10 -4 and 10 -3 , respectively. The iterative procedure terminates when the inequalities r

(k+1) pri 2 < ǫ (k+1) pri and r (k+1) dual 2 < ǫ (k+1)
dual are satisfied. Regarding the choice of the step sizes ρ 1 and ρ 2 , even though the ADMM theory ensures convergence for arbitrary positive values, it is well known that the speed of convergence will be influenced by this choice. Thus, we found a crude heuristics (ρ 1 = Ω 1 2 and ρ 2 = √ λ Ω 2 2 ) under which the computational time is well-balanced for the wide range of weights λ. Note that the classical matrix norm • 2 subordinate to the corresponding vector norm is usually estimated and then it does not involve heavy computations.

In addition, one can argue that the most expensive step of CoRE algorithm is solving the unconstrained problem [START_REF] Mohr | Comparing iterative solvers for linear systems associated with the finite difference discretisation of the forward problem in electroencephalographic source analysis[END_REF]. However, exact solution of this subproblem is not necessary to obtain the final estimate, hence this step can be efficiently executed by several runs of some appropriate iterative solver (e.g. conjugate gradient). Moreover, these can usually be warm-started, which is exactly the way we implemented the minimization (using v Null computed in the previous iteration as a starting point). It is also important to note that all matrix operators involved are extremely sparse with respect to their dimensions (the numerical density grows only linearly with the dimensions), which is an important advantage of the cosparse approach applied in this context [START_REF] Kitic | Review of cosparse signal recovery methods applied to sound source localization[END_REF].

SIMULATIONS

In this section, the numerical performance of the CoRE algorithm is studied on simulated epileptic interictal data and compared with that of the MCE [START_REF] Uutela | Visualization of magnetoencephalographic data using minimum current estimates[END_REF], RapMUSIC [START_REF] Mosher | Source localization using Recursively Applied and Projected (RAP) music[END_REF] and 4-D-MUSIC techniques [START_REF] Albera | Brain source localization using a fourth-order deflation scheme[END_REF].

Experiments and performance criterion

Two scenarios were considered for this comparison of performance. The first one aimed at studying the influence of the SNR for N = 128 electrodes while the second one allowed us to analyse the influence of the number of electrodes for an SNR value of -6 dB.

For both scenarios, P = 3 synchronous epileptic dipoles were arranged in S: P = 3 monopoles were placed at , respectively (locations are given in centimeters). Note that the origin (O) of the head model was defined as the intersection of the O-Cz axis (z-axis), the O-T4 axis (x-axis) and the O-Fpz axis (y-axis). A physiologically-relevant model [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns?[END_REF] was used to generate the epileptic activity. It is noteworthy that this activity was the same for the three epileptic dipoles, leading to synchronous epileptic sources. On the other hand, the background activity, i.e. the activity of non-epileptic dipoles of S, was generated as Gaussian and as temporally and spatially white. Its power was controlled by a multiplicative coefficient in order to get different SNR values.

θ - 1 = [
As far as the head model is concerned, we used three nested concentric spheres with radius (cm) equal to 7, 8, 9.2 and conductivities (siemens/cm) equal to 1, 0.0667, 1, 10 -10 . One hundred twenty eight electrodes were placed on the scalp sphere using the 10-5 system [START_REF] Oostenvelda | The five percent electrode system for high-resolution EEG and ERP measurements[END_REF]. The 10-5 system1 is a logical extension of the 10-20 and 10-10 systems, enabling the use of up to 345 electrodes locations. In addition, in order to apply the FDM and compute the analysis operator A, we created a cubic grid with a 0.92 cm spacing (L 3 = 9261). Consequently, the size of A and the number Q of dipoles of S were (4168 × 4168) and 626, respectively. The quality of the source localization was quantified for each method by means of the average Root Mean Square Error (RMSE), which is defined by:

RMSE = 1 P R P p=1 R r=1 min 1≤j≤P θ (-) p -θ j (r) (16)
where R is the number of realizations fixed to 100, where θ p (r) = [θ (-) p (r), θ (+) p (r)] is the position of the p-th epileptic dipole and where θ j (r) = [ θ (-) j (r), θ (+) j (r)] is the j-th dipole location estimated during the r-th Monte Carlo realization. It is noteworthy that from one realization to another, the temporal dynamics of the Q dipoles of S were changed while the location of the three epileptic dipoles stayed unchanged.

Computer results

Figure 1 shows the RMSE criterion at the output of the four algorithms as a function of the SNR. It appears that the CoRE method is more robust with respect to the presence of noise than MCE. CoRE succeeds in localizing perfectly the three epileptic dipoles beyond -17.5 dB while the synthesis-based method does not manage to do it below 0 dB. In such a context of synchronous sources, both MUSIC-like approaches do not succeed in localizing the three epileptic dipoles whatever the SNR value is.

Figure 2 displays the RMSE criterion at the output of the four algorithms as a function of the number N of electrodes for an SNR value of -6 dB. CoRE requires at least 64 electrodes to achieve a perfect result for such an SNR value while the other algorithms do not manage to perfectly localize the three synchronous epileptic dipoles whatever the number of electrodes is.

CONCLUSION

We proposed a brain source localization method, named CoRE, which uses a physics-driven structured cosparse representation of EEG signals. Our computer results on epileptic interictal data generated using a spherical head model and a physiologically-relevant source model showed the feasability of CoRE approach and its ability to deal with synchronous dipole activities, overcoming RapMUSIC and 4-D-MUSIC in such a context. We also illustrated the better behavior of CoRE for low SNR values compared to the synthesis-based MCE technique. Forthcoming work will include an evaluation of our method on data simulated using a realistic head model but also on real data. Moreover, we will test other optimization strategies to take into account the physiological constraints on current sources.

Fig. 1 .

 1 Fig. 1. Behavior of CoRE, MCE, RapMUSIC and 4-D-MUSIC as a function of the SNR for P = 3 epileptic dipoles and N = 128 electrodes.

Fig. 2 .

 2 Fig. 2. Behavior of CoRE, MCE, RapMUSIC and 4-D-MUSIC as a function of the number of electrodes for P = 3 epileptic dipole and an SNR value of -6 dB.

  1.84, 0.92, 6.44] T , θ - 2 = [-1.84, -1.84, 6.44] T and θ -

	3 ated monopoles were placed at θ + = [1.84, 5.52, 3.68] T , respectively, while the associ-1 = [0.92, 0.92, 4.6] T , θ + 2 = [-1.84, -1.84, 4.6] T and θ + 3 = [1.84, 4.6, 2.76]

It is called the 5% system or the 10-5 system because it uses proportional distances of 5% of the total length along contours between skull landmarks, compared to the

and 10% distances used in the the 10-20 and 10-10 systems, respectively.

Appendix

Using Witwer's FDM [START_REF] Witwer | The effect of media inhomogeneities upon intracranial electrical fields[END_REF], the total current flow value, I

[k] ℓ1,ℓ2,ℓ3 , at voxel (ℓ 1 , ℓ 2 , ℓ 3 ) and at time kT e , where T e the sample period, is given by: [START_REF] Oostenvelda | The five percent electrode system for high-resolution EEG and ERP measurements[END_REF] with:

ℓ1,ℓ2,ℓ3 ) as a fourth order array, A = (A (m1,m2,m3) ℓ1,ℓ2,ℓ3

) and B = (B ℓ1,ℓ2,ℓ3 ) as third order arrays, and the Hadamard product between multi-way arrays, it is easy to implement the left hand side of equation ( 17) using matrix programming languages such as Matlab.


J L Mari

B Yven

3D high resolution seismic model with depth: A relevant guide for Andra deep geological repository project

Keywords: Permeability index, density, depth conversion, elastic impedance, elastic parameters, geological model, seismic data

In the context of deep geological disposal of high level radioactive wastes, the French National Radioactive Waste Management Agency (Andra) has conducted an extensive characterization of the Callovo-Oxfordian argillaceous rock and surrounding formations in the Eastern Paris Basin. For such characterization, a detailed 3D geological model is needed. This paper shows the procedure used for building the 3D model with depth by combining time-to-depth conversion of seismic horizons, consistent seismic velocity model and elastic impedance in time. It also shows how the 3D model is used for mechanical and hydrogeological studies. The 3D field data example illustrates the potential of the proposed depth conversion procedure for estimating a density model with depth. The geological model shows good agreement with logging data obtained at a reference well, a closest borehole located in the vicinity of the 3D seismic area. Modeling of mechanical parameters such as shear modulus, young modulus, bulk modulus indicates weak variability of these parameters which confirm the homogeneity of the Callovo-Oxfordian clay. 3D modeling of a permeability index (Ik-Seis) computed from seismic attributes (instantaneous frequency, envelope, elastic impedance) and validated at the reference borehole shows promising potential for supporting hydrogeological simulation.

Introduction

3D seismic blocks and logging data, mainly acoustic and density logs, are often used for geological model building with time data. The geological model must be then converted from time data to depth data. A geostatistical approach for time-to-depth conversion of seismic horizons is often used in many geomodeling projects. The time-to-depth conversion of the selected seismic horizons by Bayesian kriging [START_REF] Sandjivy | Depth conversion and associated uncertainties using consistent velocity model: a probabilistic unified model based on Bayesian approach[END_REF]) can be used to compute a time-to-depth conversion model at the time sampling rate (1 ms). The 3D depth conversion model and the 3D impedance block allow the computation of interval velocity and density blocks. It is also possible to extract specific attributes from seismic sections to compute a seismic index (Ik-Seis factor) leading to a better understanding of the distribution of the porous and permeable bodies [START_REF] Mari | Detection of porous and permeable formations: from laboratory measurements to seismic measurements[END_REF]. The Ik-Seis, density and velocity sections can be fruitfully combined to do a mechanical and hydrogeological study with depth.

The proposed procedure has been applied both on 2D and 3D data sets. The seismic data come from High Resolution seismic surveys conducted in the vicinity of the French National Radioactive Waste Management Agency (Andra) Meuse/Haute-Marne Underground Research Laboratory (URL) located in the Meuse department in Eastern France. . Within its scientific and experimental programs devoted to investigate the disposal of high level radioactive wastes in a Callovo-Oxfordian argillaceous rock (Cox) of the eastern Paris Basin, Andra has conducted an extensive characterization of the Cox and the Oxfordian and Dogger limestones directly located above and below the target formation. More than 50 deep wells were dedicated to the geological, hydrogeological and geomechanical characterization of the Jurassic and Triassic geological series. In addition, 2D and 3D seismic data have been recorded to build a 3D structural and stratigraphic model and to estimate the variability of the hydrogeological and geomechanical parameters of the geological formations. The main purpose of this study is to describe the results of the 3D modeling of the geometry and properties of the Cox host rock in which a repository may be built. These results will serve as the basis for further simulations to assess the repository.

After a short review of the geological setting, we describe the seismic spreads which has been used to record high resolution seismic data, the processing sequence applied to the seismic data set and the procedure used to perform the time to depth conversion of the time migrated seismic sections. Finally, we show how the density seismic sections are used to detect shaly layers and how the Ik-Seis attribute can be a qualitative indicator used to detect the presence of water-bearing layers. The results obtained are compared with logging data recorded in wells EST431 and EST433 located on the 2D seismic line 07EST10.

Geological setting

The region where Andra URL is located forms a geologically simple domain of the Paris basin, with a succession of marine limestone, marl and clay rock layers. The dip of the layers is low, around 1° to 1.5° towards the north-west. The host formation, a clay rock 155 million years old, at least 130 meters thick and located between 400 and 600 meters below ground, is referred to as the Callovo-Oxfordian argillites (Cox). The top of the formation has a greater carbonate content and includes interbedded clayey layers and carbonate rock. The Cox is more homogeneous in its central part with a clay-mineral concentration of 45-50%, which corresponds to a maximum of flooding within the depositional environment of the area. A detailed study of the spatial variability of Cox geological and physical properties is found in [START_REF] Garcia | Multivariate and spatial statistical analysis of Callovo-Oxfordian physical properties from lab and borehole logs data: Towards a characterization of lateral and vertical spatial trends in the Meuse/Haute-Marne Transposition Zone[END_REF]. This formation is embedded by two thick limestone formations of Dogger and Oxfordian age. The Dogger (Bathonian, Bajocian) corresponds to the development of carbonate platforms in the Paris Basin [START_REF] Gaumet | Fondements géologiques pour la modélisation stratigraphique des systèmes carbonates. Le jurassique moyen de l'Angleterre à la Méditerranée[END_REF][START_REF] Purser | Sédimentation et diagenèse des carbonates néritiques récents[END_REF]. The clayey formation is separated from the Dogger limestone by a regional paleo-erosion surface. The Bathonian limestone formation is the first aquifer encountered in the Jurassic series.

Above the Cox, an Oxfordian (Middle and Late) carbonate platform developed with reefs and oobioclastic facies. The Oxfordian limestone formation represents the second Jurassic aquifer. In the area investigated by seismic surveys, the outcropping formations are Kimmeridgian marls and Kimmeridgian and Tithonian limestones with a thin deposit of Cretaceous.

The readers should refer to Brigaud et al. (this issue) for a detailed description of the stratigraphic architecture of the Jurassic carbonate platforms. Figure 1 shows the geological description of the Jurassic and Triassic series at well EST433 and the location of the interpreted seismic horizons.

The structural deformations have remained small for the past 150 million years, as in the rest of the Paris basin [START_REF] Guillocheau | Meso-Cenozoic geodynamic evolution of the Paris Basin: a 3D stratigraphic constraints[END_REF][START_REF] Mari | Contribution of seismic and acoustic methods to reservoir model building[END_REF]. They are essentially limited to the Gondrecourt and Marne grabens, on the boundary of the sector studied [START_REF] André | Oxfordian sedimentary dykes: Tectonic and diagenetic implications for the eastern Paris basin[END_REF][START_REF] Rocher | Intraplate paleostresses reconstructed with calcite twinning and faulting: improved method and application to the eastern Paris Basin (Lorraine, France)[END_REF].

The geological and geophysical studies have shown that the Callovo-Oxfordian layer is regular and practically flat between these faults. Available data also confirm that the region has very low seismicity. The same processing sequence has been applied both on the 2D and 3D data sets. It is a conventional seismic sequence which includes:

Seismic acquisition and processing

• data editing,

• minimum phase conversion,

• amplitude recovery (using T 2 law),

• surface consistent amplitude compensation (source and receiver),

• de-noise and wave separation on shots,

• statics application (data shifted to floating datum),

• surface consistent deconvolution,

• velocity analysis (on 480 m grid),

• surface consistent residual statics,

• second pass velocity analysis (on 240 m grid),

• surface consistent amplitude compensation (source, receiver and offset),

• interpolation and regularization in offset planes and noise attenuation,

• velocity model updating (residual move out on 240m grid),

• QC : 60 fold CMP stack (0-1400 m offsets ) with static to the final datum (450 m MSL),

• pre-stack time migration (time shifted to the final datum: 450 m MSL),

• Q compensation,

• Noise attenuation and phase conversion (statistical to zero phase),

• Acquisition foot print attenuation and band pass filter (15,160 Hz).

Andra supplied a geological model (distribution of velocities in depth) to compute static corrections.

An up-hole survey (20 wells) was conducted to calibrate the model of static corrections.

Figure 3 shows the 2 selected seismic lines (2D line 07EST10 and the in line 405) after the processing sequence just described. For all the seismic sections and logs, the depth reference is the seismic datum plane (450 m MSL). The selected lines will be used to illustrate the different steps of the procedure and the results obtained.

During the migration process [START_REF] Robein | Velocities, Time-imaging and Depth-imaging in Reflection Seismics[END_REF], the data are migrated and sorted in common image gathers (CIG) to update the velocity model before stacking. Each image gather is composed of migrated seismic traces which are functions of time (time migration) and offset or angle. The offsetangle conversion can be done during the migration process or after by using the velocity model. The common image gather after offset-angle conversion are used to perform Amplitude versus Angle analysis [START_REF] Castagna | Offset-dependent reflectivity -Theory and practice of AVO anomalies[END_REF][START_REF] Walden | Making AVO sections more robust[END_REF] and elastic inversion (EI). The amplitudes of small angle (near-offset or intercept) migrated stacks relate to change in acoustic impedance (AI) and can be inverted back to AI using post-stack inversion algorithm (acoustic inversion). The acoustic impedance is a simple function of P-wave velocity Vp and density ρ (AI = ρ.Vp). The elastic impedance (EI) is a function of P-wave velocity Vp, S-wave velocity Vs, density ρ and incidence angle [START_REF] Connolly | Elastic impedance[END_REF].

The amplitudes of angle migrated stacks can be inverted back to impedances Ip (Ip = AI = ρ.Vp) and

Is (Is = ρ.Vs) using a linearization of Zoeppritz equations for P-wave reflectivity [START_REF] Shuey | A simplification of Zoeppritz equations[END_REF]. Such a processing refers to elastic inversion [START_REF] Whitcombe | Extended elastic impedance for fluid and lithology prediction[END_REF]. A model-based elastic inversion (a priori impedance model obtained from well data), applied to the angle migrated stacks, provides 3D impedance blocks (Ip and Is blocks). The acoustic impedance sections (Ip-sections) are shown in Figure 4 for the selected seismic profiles.

A geostatistical approach for time-to-depth conversion of seismic horizons is often used in many geo-modelling projects. From a geostatistical point of view, the time-to-depth conversion of seismic horizons is a classical estimation problem involving one or more secondary variables.The converted depth and associated uncertainty can be estimated using a kriging method which can be constrained by the well markers, velocity model and interpreted horizons. For the multilayer case, the kriging estimator should take into account all the relationships between horizons determined by the velocity model associated to each layer.

The more appropriate kriging method for this problem is the Bayesian Kriging (BK) [START_REF] Sandjivy | Depth conversion and associated uncertainties using consistent velocity model: a probabilistic unified model based on Bayesian approach[END_REF][START_REF] Omre | The Bayesian bridge between simple and universal kriging[END_REF][START_REF] Omre | Bayesian kriging merging observations and qualified guesses in kriging[END_REF]. Bayesian approach provides an excellent estimator which is more general than the traditional kriging with external drift(s) and fits very well to the The input information for BK are:

• Two-way-time (TWT) maps for interpreted horizons,

• Well markers for each horizon,

• Prior velocity model and associated uncertainty for each layer,

• As any Kriging based estimator, the BK provides:

• The estimated variable (estimated depth for each horizon),

• Variance of estimation (associated uncertainty of estimated depth).

The use of BK in depth conversion has the advantage to combine the prior knowledge of the velocity model with a certain degree of uncertainty and the well data. All sources of uncertainty (velocity and time) are integrated in a consistent way in a unique probabilistic model used for estimation or simulation.

For each selected horizon, the BK provides its estimated depth Z associated with its time t. The " Z versus t " data set is interpolated in the whole space (3D block) at the time sampling rate (1 ms) in order to obtain a time to depth conversion model, using the impedance sections to estimate the short wavelength variations of the velocity model.

The time to depth conversion procedure is illustrated via the line 07EST10 (see location map, Figure 2). 10 seismic horizons numbered from 1 to 10 have been picked in time and depth converted. The 10 seismic horizons are:

1. Top of Kimmeridgian White Limestone Figure 5 shows the selected times of the 10 seismic horizons, the time-to-depth conversion model, the depth conversion of the 10 horizons and the results of the depth conversion for the seismic trace located at the CMP 654 located near the well EST433. Figure 6 shows the two amplitude seismic sections in depth.

The time-to-depth conversion model can be used to convert in depth any type of seismic sections (amplitude, velocity, density, shear modulus, Ik-Seis factor ...).

Structural study of the Callovo-Oxfordian formation

Geomechanical study

The 3D depth conversion model is used for computing of an interval P-wave velocity block (at a thin sampling rate 1 m) which is compared with the acoustic impedance block to estimate a density block. Figure 13 (top) shows the following results:

• Between the top of Lower Oxfordian and the top of carbonated Dogger: It is the depth interval of the Cox formation. On the histogram, one single mode is observed. The average density is 2.34 g/cm 3 and the associated standard deviation is 0.03 g/cm 3 .

• Between the top and the base of carbonated Dogger: One single mode is observed on the histogram. The average density is 2.56 g/cm 3 and the associated standard deviation is 0.04 g/cm 3 .

• Between the top of White Kimmeridgian limestone and the top of Beaumont dolomites (over a depth thickness of 1200 m), three modes can be observed on the histogram.

In the Callovo-Oxfordian, the highest value of the density is 2.38 g/cm 3 . A threshold of 2.38 g/cm 3 has been applied to the density section to compute a shale indicator which equals 0 if the density value is smaller than 2.38 g/cm 3 and 1 if not. The shale indicator sections are shown in Figure 13 (bottom).

The shale indicator obtained at CMP 654 has been compared with a good agreement with the indicator in the lithological log of well EST433 (Figure 1). This result allows using the shale indicator sections to study the spatial lithological variability of the shale layers from the Triassic to the Upper Jurassic formations.

Hydrogeological study

The impedance model (Ip block) can be converted into porosity by using an empirical relationship between porosity and acoustic impedance established at well locations. To model porosities, another option is to use porosity at wells location and interpolate between the wells by means of kriging. Partly due to the small number of wells, this outcome is really smooth and usually does seem geologically consistent. More dense information can be integrated in order to improve the estimation of porosity. As porosity is linked to acoustic impedance, it is relevant to use dense seismic acoustic impedance information. So a collocated cokriging of porosity integrating the seismic information is performed by using the normalized acoustic impedance as the secondary variable [START_REF] Bourges | A practical review of geostatistical processing applied to geophysical data: methods and applications[END_REF]. 3D cube makes it possible to provide 3D imaging of the connectivity of the porous bodies (Mari and Delay, 2011). Core analysis is usually carried out to establish porosity vs. permeability laws [START_REF] Zinszner | A geoscientist's guide to petrophysics[END_REF]. It has been shown that it is possible to extract new attributes from seismic sections, leading to a better understanding of the distribution of the porous and permeable bodies [START_REF] Mari | Detection of porous and permeable formations: from laboratory measurements to seismic measurements[END_REF]. The attributes are also used to detect the impermeable layers.

Laboratory experiments [START_REF] Morlier | Atténuation des ondes élastiques dans les roches poreuses saturées[END_REF] have shown that the acoustic attenuation of a clean formation can be expressed in terms of three structural parameters: porosity, permeability and specific surface. Both theoretical and experimental studies have identified the relation between acoustic attenuation and petrophysical parameters:

( )( ) 3 1 . . . 2 . . µ ρ π ϕ δ f f k S C = (1)
With δ attenuation (dB/cm), f frequency (Hz), ρf fluid density, μ fluid viscosity (centipoise) φ porosity, S Specific surface (cm 2 /cm 3 ), C calibration coefficient and k permeability (mD). [START_REF] Fabricius | Estimating permeability of carbonate rocks from porosity and Vp/Vs[END_REF] have found that the specific surface with respect to grain volume (Sg) is apparently independent from porosity. In an attempt to remove the porosity effect on Vp/Vs and mimic a reflected φ vs log (Sg) trend, they propose to use the following relationship between porosity φ, Vp/Vs and Sg:

( ) ( ) c Vs Vp b a m Sg + + = . . . log ϕ with ( ) ϕ - = 1 S Sg (2)
where it should be observed that Sg is multiplied by m to make Sg dimensionless. To establish eq 2, Fabricius at al. ( 2007) have looked at ultrasonic data, porosity, and permeability of 114 carbonate core plugs.

In practice, the parameter Ik-Seis (Indicator (I) of permeability (k) from acoustic or seismic (Seis) data) computed from equation 1 is proportional to the permeability k.

Ik-Seis = (φ.δ/S) 3 /f = (φ/SQ) 3 /f (3) with f P-wave frequency, Q quality factor, δ attenuation, S specific surface and φ porosity.

The formation permeability indicator Ik-Seis can be obtained via the computation of 4 input data: Pwave frequency and attenuation, porosity and specific surface. The procedure has been firstly conducted in acoustic logging to estimate permeability of porous layers and to detect water inflows [START_REF] Mari | Detection of porous and permeable formations: from laboratory measurements to seismic measurements[END_REF]. In seismic, the processing is performed to measure these parameters. The analytic signal is used to compute the instantaneous frequency and attenuation (Q factor). The porosity and specific surface are computed from seismic impedances obtained by elastic inversion of the migrated seismic sections. In the domain of seismic frequencies, the Ik-Seis factor can only be seen as a relative indicator which varies from 0 for less porous and permeable bodies to 1 for more porous and permeable bodies. More information concerning the data processing and analysis is given in [START_REF] Mari | Detection of porous and permeable formations: from laboratory measurements to seismic measurements[END_REF].

Well The analysis of the acoustic waves recorded simultaneously on both receivers of the acoustic tool is used to compute additional logs defined as acoustic attributes useful for the characterization of the formation, such as amplitude, shape index and attenuation logs. The results obtained are optimum if the studied wave is extracted from the records and if the signal to noise ratio is high. We show the benefit of using Singular Value Decomposition (SVD) for that purpose [START_REF] Glangeaud | Signal Processing in Geosciences[END_REF].

The SVD processing is done on the 2 constant offset sections independently, in a 5 traces (N=5) depth running window. After flattening of each constant offset section with the picked times of the refracted wave, the refracted wave signal space is given by the first eigensection obtained by SVD:

r sig = λ 1 u 1 v 1 T (4)
v 1 is the first singular vector giving the time dependence, hence named normalized wavelet, u 1 is the first singular vector giving the amplitude in depth, therefore called propagation vector and λ 1 the associated eigenvalue. The amplitude variation of the refracted wavelet over the depth interval is λ 1 u 1. It is therefore suggested that changes in phase or distortion of the acoustic signal is linked to propagation through a porous and permeable zone. The distortions can be measured by a shape index attribute. To measure the shape variation, an acoustic attribute, named Ic, independent of the energy of the source, has been introduced (Lebreton and Morlier, 1986). The Ic parameter is given by the following equation:

Ic = ((A 2 + A 3 )/ A 1 ) n (5)
where A 1 , A 2 and A 3 are the amplitudes of the first three arches, respectively, of the studied signal and n an exponent.

The shape index is computed from formula 5 with an exponent value of 3 (n = 3). shown in Figure 16 (bottom). The shaly layers, associated with a shale index equal to 0 (Figure 13) have the lowest Ik-Seis value. The Ik-Seis sections show that the Callovo-Oxfordian is a very low permeable formation with a low variability of the Ik-Seis factor (average value: 0.06, standard deviation: 0.1). Porous and permeable bodies in the Oxfordian have been detected and confirmed by acoustic logging. In the two Ik-Seis seismic sections, the Oxfordian limestone appears to have a higher Ik-Seis factor than the Dogger limestone. In the Dogger, the Ik-Seis factor is generally low except directly above the base of the target interval. However this level of higher Ik-Seis seems to vary spatially. Figure 17 is a comparative display of Ik-Seis factor at CMP 654, free fluid NMR porosity, gamma-ray and sonic logs at the well EST433. A good correlation between the Ik-Seis factor, the NMR porosity and the location of the porous levels (Hp1-Hp4 in the Oxfordian, HpD1 and HpD2 in the Dogger) can be noted. This result allows to perform a detailed study the 3D spatial variability of the Ik-Seis factor in the two carbonate formation.

Conclusions

We have described the workflow which was used to build a geological model with depth combining time-to-depth conversion of seismic horizons, consistent seismic velocity model and elastic impedance inversion in time. The field data example illustrates the potential of the proposed depth conversion procedure for estimating a density model in depth. The knowledge of the density distribution allows the computation of mechanical parameters such as shear modulus, Young modulus or bulk modulus. The seismic results obtained in the vicinity of well EST433 have been compared with the information obtained by logging (mainly velocity and density measurements). The correlation coefficients between the seismic and logging data are high (larger than 0.9) for the velocities. If the results (time depth conversion, velocities, density, mechanical parameters) remain in the expected range of uncertainties, the proposed workflow could be used to build a reliable 3D model with depth from the full seismic 3D block.

On the studied seismic lines (07EST10 and IL-405), the variability of the seismic parameters (Pwave and S-wave velocities, density, shear modulus, Young modulus and bulk modulus) is weak in the Callovo-Oxfordian leading us to conclude that it is an homogeneous formation.

The Ik-Seis sections in depth have shown that the Callovo-Oxfordian is a very low permeable formation and have detected porous and permeable bodies in the Oxfordian and at the top of the Dogger limestone. The results have been confirmed at well (EST431) using full waveform acoustic data.

A permeability index in 3D (Ik-Seis block) has a promising potential for supporting hydrogeological assessment such as determining and mapping in 3D the porosity and permeability variability of certain geological formations and hydrogeological modeling.

Figure 2

 2 Figure 2 shows the 3D seismic data acquisition area and the location of the selected seismic lines in the transposition area: the 2D line 07EST10 and the in line 405 extracted from the 3D block. The 2D line 07EST10 is oriented W-E and the in line 405 is oriented NNW-SSE. No wells being located in the 3D area, three additional 2D lines have been recorded to calibrate the 3D data on 3 wells situated outside of the 3D area. Well EST433 is located on one of the additional 2D line 07EST10 in the vicinity of the CMP 654 (Common Mid Point). The 3D design is a cross spread. The active spread is composed of 12 receiver lines with 120 stations each. The source lines are perpendicular to the receiver lines. The receiver and source line spacing's are respectively 80 m and 120 m. The receiver and source point spacing's are 20 m. The source is a vibroseis source generating a signal in the 14 -140 Hz frequency bandwidth. The bin size is 10 x 10 m. The nominal fold is 60. The in line 405 is composed of 544 CMP points. The 2D design is a split dip spread, composed of 120 stations. The receiver and source point spacing's are 25 m. The source is a vibroseis source generating a signal in the 14 -140 Hz frequency bandwidth. The bin size is 12.5 m. The nominal fold is 60. The selected part of the 2D line 07EST10 is composed of 727 CMP points.

  needs for time-to-depth conversion of seismic horizons. The advantage of BK compared to other estimators consists in the fact that the uncertainty on the trend velocity model and the local uncertainty defined by the uncertainty of interpreted time maps and local fluctuations of interval velocities can be managed simultaneously.

  Local uncertainty definition for each time map (picking uncertainty, and spatial variogram definition), • Local uncertainty definition of interval velocity for each layer (local velocity fluctuations around the trend model, and spatial variogram definition).

  2. Top of Porous Horizon HP4 3. Top of Lower Oxfordian (Top of target interval) 4. Top of Upper Callovian (RIO) 5. Top of Carbonated Dogger (Base of target interval) 6. Base of Argillaceous limestone and marls 7. Base of Carbonated Dogger 8. Top of Domerian 9. Base of Lias (base of Gryphees limestone) 10. Top of Beaumont dolomite.

Figure 7

 7 Figure7shows the elevation map below sea level of the base of the Cox formation obtained directly from the seismic data interpretation and the time-to-depth conversion process (horizon 5). The isolines are gently dipping towards the north-west (center of Paris Basin) and are subparallel. This dip has been computed and shows an average value of 1°. The difference of elevation from the east to the west in the ZIRA is 110 m. No evidence of fault may be interpreted from this map. The same conclusions arise from the maps directly below and above the Cox formation.

Figure 8

 8 Figure8shows the thickness of the Cox formation. The thickness increases from the south-west (142 m) to the north-east (162 m) in the ZIRA with a mean thickness of 153 m. The direction of thickening is coherent with the main source and direction of sediments transport.

Figure 9

 9 Figure9shows the P-wave velocity and density distributions in depth. The low values of density and velocity in the 600-750 m depth interval for the line IL-405 and in the 500-700 m depth interval for the line 07EST10 are associated with the target interval situated between the top of Lower Oxfordian and the top of the Carbonated Dogger. The knowledge of the density (ρ) and velocity (Vp, Vs) distributions allows the computation of mechanical modules such as shear modulus, Young modulus, bulk modulus and Poisson's ratio. Figure10shows the S-wave velocity and shear modulus distributions in depth.The seismic data obtained at CMP 654 have been compared with logging data recorded at well EST433.

  EST 433 was used to study formations ranging from Oxfordian to Trias. Well EST 431 was drilled to complement the geological and hydrogeological knowledge of the Oxfordian formation. This formation consists essentially of limestone deposited in a vast sedimentary platform. The limestone facies, which vary from one borehole to another, are generally bio-detritic with reef constructions. In this formation, porosity ranges between 5 and 20% and "porous horizons" of kilometric extension have been identified. The observed water inflows are usually located in high porosity zones. During the drilling, water inflows were detected at -456 m and -503 m (reference depth: seismic datum plane).The acoustic tool used for the field experiment is a flexible monopole tool with two pairs of receivers: a pair of near receivers (1 and 1.25 m offsets) and a pair of far receivers (3 m (R1) and 3.25 m (R2) offsets). The data have been recorded through the far offset configuration. The sampling depth interval is 10 cm. The sampling time interval is 5 microseconds. The length of recording is 5 ms. The acoustic log has been run in the Oxfordian carbonate formation, in the 420-600 m depth interval. Figure14(top left) shows the R1 and R2 constant offset sections in the 0.5-1.25 ms time interval where the refracted P-waves appear. The processing of the acoustic data has been described in detail byMari et al. (2011). Figure14(top right and bottom left) is a display of acoustic logs: P-wave velocity (V P ), P-wave frequency, acoustic porosity, P-wave attenuation. For a clean formation, if the matrix and fluid velocities are known, an acoustic porosity log can be computed from the acoustic V p velocities using the formula given by[START_REF] Wyllie | Elastic wave velocities in heterogeneous and porous media[END_REF] expressed in velocities. The matrix velocity value has been chosen at 6300 m/s, and the fluid velocity at 1500 m/s (the formation fluid is water). The acoustic porosity log, valid only in the clean part of the formation, shows a strong correlation (correlation coefficient: 0.86) with a NMR porosity log (not displayed here) recorded in the well. Figure14(bottom right) shows a comparison between acoustic velocities at wells EST431 (black curve) and EST433 (red curve) and between acoustic velocities at well EST431 (black curve) and seismic velocities at CMP 654 (red curve). The correlation coefficients are high (0.88). The attenuation (expressed in dB/m) of the formation is computed from the first eigensection (obtained by Singular Value Decomposition SVD) of the refracted P-wave acoustic signal recorded by the two adjacent receivers of the acoustic tool. The results obtained by the SVD processing procedure are shown in Figure 15 (top). The Ik-Seis log detects three permeable zones at 456 m, between 490 and 530 m, and at 594 m. The permeable zone located at 594 m corresponds to a high value of conductivity and is characterized by a low porosity (6 %), a 10 dB/m attenuation, but a significant decrease of the P-wave frequency. The hydraulic tests and conductivity measurements conducted later on did not confirm the inflow at 456 m seen during the drilling, but have validated the 490-530 m and 594 m permeable zones detected by the acoustic logging.

Figure 15 (

 15 Figure 15 (top) shows the normalized wavelet (v 1 ) and the associated amplitude (λ 1 u 1 ) log versus depth, for the two constant offset sections associated with the two receivers (R1 and R2) of the acoustic tool. The amplitude logs have been used to compute the attenuation log (Figure 14, bottom left) expressed in dB/m. The correlation coefficient (Figure 15, bottom left) between the two normalized wavelets has been computed at each depth. We can notice some anomalies at local depth (446, 465, 550, 580 and 590 m) and a significant decrease of the correlation coefficient in the 490-530 m depth interval. The interval corresponds to the porous and permeable zone detected by the Ik-Seis factor (Figure 15, bottom right). It is therefore suggested that changes in phase or distortion of the acoustic

  Figure 15 (top) allows the comparison between the two shape index logs (Ic-R1 and Ic-R2). The shape index logs highlight anomalic zones, in the 440 -470 m depth interval with a maximum value at 454 m, and in the 490 -530 m depth interval with a peak at 495 m which corresponds to the porous layer HP4. In order to reduce the noise to extract the common component of the two shape index logs (Ic-R1 and Ic-R2), the geometric mean of the two has been computed. The resulting shape index log Ic is compared with the Ik-Seis log and the correlation coefficient log (Figure15, bottom left). The comparison shows a good coherence between the Ic log and the correlation coefficient log. The permeable and porous layers in the 490 -530 m depth interval are seen both by the shape index log and by the Ik-Seis log.

Figure 15 (

 15 Figure 15 (bottom left) shows a comparison between Ik-Seis factors computed from acoustic data (well EST431, red curve) and seismic data (CMP 654, black curve). The seismic Ik-Seis factor is the envelope of the acoustic Ik-Seis factor, due to the different of vertical resolution between the seismic and acoustic data. The seismic Ik-Seis factor points out to permeable zones in the 440 -470 m and 490 -530 m depth intervals. The flow at 594 m cannot be detected due to the lack of resolution. The vertical resolution is proportional to the seismic wavelength; a quarter of the wavelength is usually adopted. The seismic velocity blocks (Figure 9) and the instantaneous frequencies have been used to compute the vertical resolution seismic sections (Figure 16, top). The Ik-seis seismic sections are

Figure 1 :

 1 Figure 1: EST 433 geological log and seismic horizons. 1 Top of Kimmeridgian White Limestones, 2 Top of Porous Horizon HP4, 3 Top of Lower Oxfordian (Top of target interval), 4 Top of Upper Callovian (RIO), 5 Top of Carbonated Dogger (Base of target interval), 6 Base of Argillaceous limestones and marls, 7 Base of Carbonated Dogger, 8 Top of Domerian, 9 Base of Lias (base of Gryphees limestone, 10 Top of Beaumont dolomite.

Figure 2 :

 2 Figure 2: Location map: 3D seismic area, seismic lines (IL 405 and 07EST10) and Well EST433. The 2D line 07EST10 is oriented W-E and the in line 405 is oriented NNW-SSE.

Figure 3 :

 3 Figure 3: Amplitude seismic sections in time.

Figure 4 :

 4 Figure 4: Acoustic impedance (g/cm 3 m/s) sections in time.

Figure 5 :

 5 Figure 5: Time-to-depth conversion by Bayesian Kriging.

Figure 6 :

 6 Figure 6: Amplitude seismic sections in depth.

Figure 7 :

 7 Figure 7: Elevation map of the base of the Callovo-Oxfordian formation.

Figure 8 :

 8 Figure 8: Thickness map of the Callovo-Oxfordian formation.

Figure 9 :

 9 Figure 9: P-wave velocity and density seismic sections in depth.

Figure 10 :

 10 Figure 10: S-wave velocity and shear modulus seismic sections in depth.

Figure 11 :

 11 Figure 11: Comparison between seismic and logging data at well EST433 (CMP 654, line 07EST10).

Figure 12 :

 12 Figure 12: Cross plots between seismic (CMP 654, line 07EST10) and logging data at well EST433.

Figure 13 :

 13 Figure 13: Density analysis and shale indicator sections in depth

Figure 15 :

 15 Figure 15: Permeability estimation from acoustic logs at well EST431. Top: Analysis of refracted waves by SVD (Wavelet, amplitude, shape index) on receivers R1 and R2 of the acoustic logging tool. Bottom left: predicted permeability (Ik-Seis), shape index and correlation coefficient logs. Bottom right: comparison between Ik-Seis factors computed from acoustic data (well EST431) and seismic data (CMP 654).

Figure 16 :

 16 Figure 16: Vertical resolution and Ik-Seis seismic sections in depth.

Figure 17 :

 17 Figure 17: Correlation between permeability index Ik-Seis and logging data (NMR porosity, Gammaray and sonic logs).

  Figure 1

  Figure 5

Figure 9 Figure 10 Figure 11

 91011 Figure 9

Figure 12

 12 Figure 12

Figure 13 Figure 14

 1314 Figure 13

Figure 15

 15 Figure 15

Figure 16

 16 Figure 16

Figure 17

 17 Figure 17

Table 1 :

 1 Variability of mechanical parameters in the Callovo-Oxfordian formation

		Average Value	Standard deviation
	P-wave velocity (m/s)	3000	210
	S-wave velocity (m/s)	1700	200
	Density (g/cm 3 )	2.34	0.03
	Poisson's ratio	0.28	0.06
	Bulk modulus (GPa)	13	2.7
	Shear modulus (GPa)	6.5	1.5
	Young modulus (GPa)	17	3.1

Acknowledgements

This work was financed by Andra. We thank Michel Hayet and Philippe Dubreuilh for their comments and suggestions. The authors acknowledge Arben Shtuka (Seisquare) for very useful discussions on various occasions, specifically for his experience in the use of geostatistical methods for time to depth conversion.


kernel H(x, v)=a(|x|) (-v).

8 > > < > > : ẋi = v i 2 R d vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i ) 2 R d , for i =1,...N,
where a(r): =a β (r)= 1 (1+r 2 ) β , β>0 models the exchange of information between agents.

Mattia Bongini

(Un)conditional consensus emergence under feedback controls 5o f2 8

The Cucker-Smale model is obtained by the choice of the interaction kernel H(x, v)=a(|x|) (-v).

8 > > < > > : ẋi = v i 2 R d vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i ) 2 R d , for i =1,...N,
where a(r): =a β (r)= 1 (1+r 2 ) β , β>0 models the exchange of information between agents. We say that a solution {(x i ,v i )} N i=1 tends to consensus if there exists a vector w 2 R d such that lim t!+1 kv i (t) -wk = 0 for every i =1,...,N.

Mattia Bongini

(Un)conditional consensus emergence under feedback controls 5o f2 8

Consider the symmetric bilinear form

B(u, v)= 1 2N 2 X i,j hu i -u j ,v i -v j i = 1 N N X i=1 hu i ,v i i-hū, vi,
and X(t)=B(x(t),x(t)),V (t)=B(v(t),v(t)).

Mattia Bongini (Un)conditional consensus emergence under feedback controls 6o f2 8

Conditional consensus emergence for generic a(•)

Consider the symmetric bilinear form

B(u, v)= 1 2N 2 X i,j hu i -u j ,v i -v j i = 1 N N X i=1 hu i ,v i i-hū, vi,
and

X(t)=B(x(t),x(t)),V (t)=B(v(t),v(t)).
Theorem (Ha-Ha-Kim)

Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N be such that X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ) satisfy Z 1 p X 0 a( p 2Nr)dr ≥ p V 0 .
Then the solution with initial data (x 0 ,v 0 ) tends to consensus.

Mattia Bongini (Un)conditional consensus emergence under feedback controls 6o f2 8

Conditional consensus emergence for generic a(•)

Consider the symmetric bilinear form

B(u, v)= 1 2N 2 X i,j hu i -u j ,v i -v j i = 1 N N X i=1 hu i ,v i i-hū, vi,
and

X(t)=B(x(t),x(t)),V (t)=B(v(t),v(t)).
Theorem (Ha-Ha-Kim)

Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N be such that X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ) satisfy Z 1 p X 0 a( p 2Nr)dr ≥ p V 0 .
Then the solution with initial data (x 0 ,v 0 ) tends to consensus. Hence for a(t): =a β (t) we have unconditional convergence for 0  β  1 2 .

If a(•)=a β (•) and β>1/2 then the consensus condition is not satisfied by all (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N , as shown by Cucker and Smale.

Consider β = 1, N = 2, d = 1 and let x(t)=x 1 (t) -x 2 (t) and v(t)=v 1 (t) -v 2 (t) be the relative positions and velocities of the two agents on the line: rewrite the system as 8 <

:

ẋ = v v = - v 1+x 2 with initial conditions x(0) = x 0 and v(0) = v 0 > 0.
Mattia Bongini (Un)conditional consensus emergence under feedback controls 7o f2 8

If a(•)=a β (•) and β>1/2 then the consensus condition is not satisfied by all (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N , as shown by Cucker and Smale.

Consider β = 1, N = 2, d = 1 and let x(t)=x 1 (t) -x 2 (t) and v(t)=v 1 (t) -v 2 (t) be the relative positions and velocities of the two agents on the line: rewrite the system as 8 <

:

ẋ = v v = - v 1+x 2 with initial conditions x(0) = x 0 and v(0) = v 0 > 0. By direct integration v(t)=-arctan x(t) + arctan x 0 + v 0 .
Hence, if arctan x 0 + v 0 >⇡/2+" we have v(t) >⇡/2+" -arctan x(t) >", 8t 2 R + .

Mattia Bongini

(Un)conditional consensus emergence under feedback controls 7o f2 8

However, in the work Caponigro-Fornasier-Piccoli-Trélat it is shown that the addition of an extra term of the form

8 > > < > > : ẋi = v i vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i )+γ(v -v i ), for i =1,...N,
where γ>0, implies convergence to consensus unconditionally.

Mattia Bongini (Un)conditional consensus emergence under feedback controls 8o f2 8

However, in the work Caponigro-Fornasier-Piccoli-Trélat it is shown that the addition of an extra term of the form

8 > > < > > : ẋi = v i vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i )+γ(v -v i ), for i =1,...N,
where γ>0, implies convergence to consensus unconditionally. This happens because each agent has an additional stirring to reach the common goal v.

But what happens if every agent has its own "local" goal?

Mattia Bongini (Un)conditional consensus emergence under feedback controls 8o f2 8

Introducing decentralized controls

We study how the convergence to consensus is modified by the introduction of the following extra term

8 > > < > > : ẋi = v i vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i )+γ(v R i -v i ), for i =1,...N,
where γ>0 and

v R i = 1 P N k=1 χ [0,R] (kx i -x k k) N X j=1 χ [0,R] (kx i -x j k)v j ,
is the local mean velocity calculated by the i-th agent taking into account only the agents inside a ball of radius R centered on him.

Mattia Bongini (Un)conditional consensus emergence under feedback controls 9o f2 8

Introducing decentralized controls

We study how the convergence to consensus is modified by the introduction of the following extra term

8 > > < > > : ẋi = v i vi = 1 N N X j=1 a ⇣ kx i -x j k 2 ⌘ (v j -v i )+γ(v R i -v i ), for i =1,...N,
where γ>0 and

v R i = 1 P N k=1 χ [0,R] (kx i -x k k) N X j=1 χ [0,R] (kx i -x j k)v j ,
is the local mean velocity calculated by the i-th agent taking into account only the agents inside a ball of radius R centered on him. GOAL: Study the behavior of the system for R ! 0 and R ! +1.

Mattia Bongini (Un)conditional consensus emergence under feedback controls 9o f2 8

Behavior of the modified system in function of R To avoid the issues coming from a normalizing factor varying from agent to agent, in our study of the asymptotic behavior we replace it with a function ⌘ R (t) 2 L 1 ([0, +1)) satisfying

lim R!0 ηR(t)=1 =l i m R!0 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , lim R!+1 ηR(t)=N =l i m R!+1 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , Mattia Bongini (Un)conditional consensus emergence under feedback controls
To avoid the issues coming from a normalizing factor varying from agent to agent, in our study of the asymptotic behavior we replace it with a function ⌘ R (t) 2 L 1 ([0, +1)) satisfying

lim R!0 ηR(t)=1 =l i m R!0 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , lim R!+1 ηR(t)=N =l i m R!+1 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , like ηR(t)= max i=1,...,N ( N X k=1 χ [0,R] (kxi(t) -x k (t)k)
)

.

To avoid the issues coming from a normalizing factor varying from agent to agent, in our study of the asymptotic behavior we replace it with a function ⌘ R (t) 2 L 1 ([0, +1)) satisfying

lim R!0 ηR(t)=1 =l i m R!0 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , lim R!+1 ηR(t)=N =l i m R!+1 N X k=1 χ [0,R] (kxi(t) -x k (t)k) for every i ! , like ηR(t)= max i=1,...,N ( N X k=1 χ [0,R] (kxi(t) -x k (t)k)
)

.

Hence, we replace v R i in the system with

vR i = 1 ηR N X k=1 χ [0,R] (kxi -x k k)vj. Theorem Replace v R i with vR i in the model. Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N
and define X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ). Then the solution with initial data (x 0 ,v 0 ) tends to consensus if one of the following is satisfied:

⌅ p 2NX 0  R and R 1 p X 0 a( p 2Nr)dr + γ ⇣ R p 2N - p X 0 ⌘ ≥ p V 0 ; Theorem Replace v R i with vR i in the model. Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N
and define X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ). Then the solution with initial data (x 0 ,v 0 ) tends to consensus if one of the following is satisfied:

⌅ p 2NX 0  R and R 1 p X 0 a( p 2Nr)dr + γ ⇣ R p 2N - p X 0 ⌘ ≥ p V 0 ; ⌅ p 2NX 0 >Rand R 1 p X 0 a( p 2Nr)dr ≥ p V 0 . Theorem Replace v R i with vR i in the model. Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N
and define X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ). Then the solution with initial data (x 0 ,v 0 ) tends to consensus if one of the following is satisfied:

⌅ p 2NX 0  R and R 1 p X 0 a( p 2Nr)dr + γ ⇣ R p 2N - p X 0 ⌘ ≥ p V 0 ; ⌅ p 2NX 0 >Rand R 1 p X 0 a( p 2Nr)dr ≥ p V 0 .
Proof.

via the Ha-Ha-Kim Theorem. In this case, the Lyapunov functional V decreases at least as

V (t)  V (0)e γt .
Mattia Bongini (Un)conditional consensus emergence under feedback controls ⌅ We have shown that an unconditional (w.r.t. the initial data) convergence to consensus is guaranteed only in the case R = 1, i.e., the agents are completely "open-minded". In this case, the Lyapunov functional V decreases at least as

V (t)  V (0)e γt .
⌅ Question: are there more economical and efficient strategies which guarantees unconditional stabilization?

Mattia Bongini (Un)conditional consensus emergence under feedback controls ⌅ We have shown that an unconditional (w.r.t. the initial data) convergence to consensus is guaranteed only in the case R = 1, i.e., the agents are completely "open-minded". In this case, the Lyapunov functional V decreases at least as

V (t)  V (0)e γt .
⌅ Question: are there more economical and efficient strategies which guarantees unconditional stabilization?

⌅ More precisely, we wish to make the quantity d dt V (t) the smallest possible and use the minimal amount of intervention (e.g., by imposing sparsifying constraints).

Mattia Bongini

(Un)conditional consensus emergence under feedback controls The 1-sparse feedback strategy Definition For every t ≥ 0w ed e fi n et h e1-sparse feedback control as

ui(t)= ( -M v ⊥ i (t) kv ⊥ i (t)k if i = b ι(t) 0i f i 6 = b ι(t)
where b ◆(t) 2{1,...,N} is the minimum index such that

kv ⊥ b ι(t) (t)k = max j=1,...,N kv ⊥ j (t)k. Mattia Bongini (Un)conditional consensus emergence under feedback controls > > > < > > > : ẋi = vi 2 R d vi = -bivi + N X j=1 a kxi -xjk 2 (xj -xi)+ N X j=1 j6 =i f kxi -xjk 2 (xi -xj) 2 R d
where ] is the friction acting on the system, Mattia Bongini (Un)conditional consensus emergence under feedback controls > > > < > > > :

⌅ b i :[0, +1) ! [0, Λ
ẋi = vi 2 R d vi = -bivi + N X j=1 a kxi -xjk 2 (xj -xi)+ N X j=1 j6 =i f kxi -xjk 2 (xi -xj) 2 R d
where

⌅ b i :[0, +1) ! [0, Λ] is the friction acting on the system, ⌅ a :[0, +1) ! [0, +1) is the rate of communication, Mattia Bongini (Un)conditional consensus emergence under feedback controls > > > < > > > : ẋi = vi 2 R d vi = -bivi + N X j=1 a kxi -xjk 2 (xj -xi)+ N X j=1 j6 =i f kxi -xjk 2 (xi -xj) 2 R d
where under the sparsifying constraint P N i=1 ku i (t)kM for M>0.

⌅ b i :[0, +1) ! [0, Λ] is the friction acting on the system, ⌅ a :[0, +1) ! [0, +1) is the rate of communication, ⌅ f :(0, +1) ! (0, +1) such that Z +1 δ f (

  conditional consensus emergence under feedback controls

⌅

  We have shown that an unconditional (w.r.t. the initial data) convergence to consensus is guaranteed only in the case R = 1, i.e., the agents are completely "open-minded".

For

  every initial condition (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N and M>0 there exist T>0 and a sparse control u :[ 0 ,T] ! (R d ) N ,w i t h P N i=1 ku i (t)kM for every t 2 [0,T] such that the associated solution tends to consensus.Mattia Bongini(Un)conditional consensus emergence under feedback controls

For

  every initial condition (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N and M>0 there exist T>0 and a sparse control u :[ 0 ,T] ! (R d ) N ,w i t h P N i=1 ku i (t)kM for every t 2 [0,T] such that the associated solution tends to consensus. We can choose adaptively the control law explicitly as one of the solutions of the variational problem min B(v, u)+ T (threshold function byHa-Ha-Kim)

  r) dr < 1 for every δ>0, Difference f (r)a(r) for Lennard-Jones potentials.

⌅

  Assume we are in the case β>1 and E(0) ≥ #. Can we again stabilize the society by external parsimonious intervention? ⌅ We again introduce a control term inside the model xjk 2 (xixj) + ui

Mattia Bongini(Un)conditional consensus emergence under feedback controls

a(r)= σ a r 7 and f (r)= σ f r 13 .Mattia Bongini (Un)conditional consensus emergence under feedback controls

⌅ the kinetic energy K(t):= 1 2P N i=1 kv i (t)k 2 ,Mattia Bongini (Un)conditional consensus emergence under feedback controls

⌅ Assume we are in the case β>1 and E(0) ≥ #. Can we again stabilize the society by external parsimonious intervention?Mattia Bongini (Un)conditional consensus emergence under feedback controls

⌅ the 1-sparse feedback control belongs to the set argmin u2K(t) J t,q (u)for every t  T ;Mattia Bongini (Un)conditional consensus emergence under feedback controls

⌅ We have shown a quantitative link between a degree of "narrowmindedness" and convergence to global patterns;Mattia Bongini (Un)conditional consensus emergence under feedback controls

⌅ We have shown a quantitative link between a degree of "narrowmindedness" and convergence to global patterns;

⌅ a sparse control strategy is more efficient than local decentralized control in order to enforce convergence: if we pay our attention solely to the "most stubborn" agent while leaving the other free to adjust themselves, we actually have the "best greedy" strategy;

⌅ in contrast to what happen with the Cucker-Smale model, the result for the 1-sparse control feedback in the Cucker-Dong model is conditional (it depends on the initial conditions of the system) ) we don't know if the conditions are necessary.

Mattia Bongini (Un)conditional consensus emergence under feedback controls

Theorem

Replace v R i with vR i in the model. Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N and define X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ). Then the solution with initial data (x 0 ,v 0 ) tends to consensus if one of the following is satisfied:

As u ffi c i e n tc o n d i t i o n Theorem Replace v R i with vR i in the model. Let (x 0 ,v 0 ) 2 (R d ) N ⇥ (R d ) N and define X 0 = B(x 0 ,x 0 ) and V 0 = B(v 0 ,v 0 ). Then the solution with initial data (x 0 ,v 0 ) tends to consensus if one of the following is satisfied:

Proof.

via the Ha-Ha-Kim Theorem.

⌅ Result generalizable to different forms of ⌘ R and χ [0,R] .

Mattia Bongini (Un)conditional consensus emergence under feedback controls

We introduce a control term in the Cucker-Smale model:

The admissible controls will be measurable functions u =(u 1 ,...,u N ): [0, +1) ! R N such that P N i=1 ku i (t)kM for every t>0, for a given constant M .

We introduce a control term in the Cucker-Smale model:

The admissible controls will be measurable functions u =(u 1 ,...,u N ): [0, +1) ! R N such that P N i=1 ku i (t)kM for every t>0, for a given constant M .

Our aim is then to find admissible controls steering the system to the consensus region and minimizing

in a large class of controls.

Mattia Bongini (Un)conditional consensus emergence under feedback controls

Definition

For every t ≥ 0w ed e fi n et h e1-sparse feedback control as

where b ◆(t) 2{1,...,N} is the minimum index such that

Theorem

The 1-sparse control is the minimizer of

among all the solution of the variational problem.

Definition

For every t ≥ 0w ed e fi n et h e1-sparse feedback control as

where b ◆(t) 2{1,...,N} is the minimum index such that

Theorem

The 1-sparse control is the minimizer of

among all the solution of the variational problem.

The control acts on the most "stubborn" agent. We may call it the "shepherd dog strategy"

The Cucker-Dong model: for every 1

Mattia Bongini (Un)conditional consensus emergence under feedback controls

Conditional consensus emergence

Theorem (Cucker -Dong)

Consider a population of N agents modeled by a Cucker-Dong system with a(t):=a

Then there exists a unique solution (x(t),v(t)) of the system with initial state (x(0),v(0)). Moreover if one of the two following hypotheses holds:

Conditional consensus emergence Theorem (Cucker -Dong)

Consider a population of N agents modeled by a Cucker-Dong system with a(t):=a

Then there exists a unique solution (x(t),v(t)) of the system with initial state (x(0),v(0)). Moreover if one of the two following hypotheses holds:

Conditional consensus emergence

Theorem (Cucker -Dong)

Consider a population of N agents modeled by a Cucker-Dong system with a(t):=a

Then there exists a unique solution (x(t),v(t)) of the system with initial state (x(0),v(0)). Moreover if one of the two following hypotheses holds:

Conditional consensus emergence

Theorem (Cucker -Dong)

Consider a population of N agents modeled by a Cucker-Dong system with a(t):=a

Then there exists a unique solution (x(t),v(t)) of the system with initial state (x(0),v(0)). Moreover if one of the two following hypotheses holds:

R 1 0 a(r)dr then the population is cohesive and collision-avoiding, i.e., there exist two constants B 0 and b 0 > 0 such that, for every t ≥ 0 b 0 kx i (t) -x j (t)kB 0 for all 1  i 6 = j  N.

For every t ≥ 0w ed e fi n et h e1-sparse feedback control as

where b ◆(t) 2{1,...,N} is the minimum index such that kv b ι(t) (t)k = max j=1,...,N kv j (t)k.

Mattia Bongini

(Un)conditional consensus emergence under feedback controls

Definition

For every t ≥ 0w ed e fi n et h e1-sparse feedback control as

where b ◆(t) 2{1,...,N} is the minimum index such that

Define for every t>0t h es e t

and for every t>0 and q>0 the functional J t,q :(R d ) N ! R Jt,q(u)=hv(t),ui +

) N be such that the following hold:

) N be such that the following hold:

Mattia Bongini (Un)conditional consensus emergence under feedback controls

) N be such that the following hold:

Then...

Mattia Bongini (Un)conditional consensus emergence under feedback controls ... then there exist T>0 and q>0 such that

⌅ the 1-sparse feedback control belongs to the set argmin u2K(t) J t,q (u) for every t  T ;

⌅ the 1-sparse feedback control is an instantaneous minimizer of the functional

over all possible feedback controls in argmin u2K(t) J t,q (u);

Mattia Bongini (Un)conditional consensus emergence under feedback controls ... then there exist T>0 and q>0 such that

⌅ the 1-sparse feedback control belongs to the set argmin u2K(t) J t,q (u) for every t  T ;

⌅ the 1-sparse feedback control is an instantaneous minimizer of the functional

over all possible feedback controls in argmin u2K(t) J t,q (u);

⌅ there exists a solution of the controlled Cucker-Dong system associated to a control solving the variational problem e u 2 argmin u∈K(t)

J t,q (u) for every t  T.

Mattia Bongini (Un)conditional consensus emergence under feedback controls

Summing up our results

⌅ We have shown a quantitative link between a degree of "narrowmindedness" and convergence to global patterns;

⌅ a sparse control strategy is more efficient than local decentralized control in order to enforce convergence: if we pay our attention solely to the "most stubborn" agent while leaving the other free to adjust themselves, we actually have the "best greedy" strategy;

Mattia Bongini (Un)conditional consensus emergence under feedback controls

Summing up our results

⌅ We have shown a quantitative link between a degree of "narrowmindedness" and convergence to global patterns;

⌅ a sparse control strategy is more efficient than local decentralized control in order to enforce convergence: if we pay our attention solely to the "most stubborn" agent while leaving the other free to adjust themselves, we actually have the "best greedy" strategy;