Céline Hadziioannou

Éric Larose

Adam Baig

Philippe Roux

Michel Campillo

E Larose

Improving temporal resolution in

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

It was observed that the two-point cross-correlation between long records of ambient seismic noise yields the Green Function (GF) between the two stations (e.g. [START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]). This property is expected when noise sources are distributed isotropically around the stations which is 1 eventually the case when the ambient noise correlation is integrated over a long time [START_REF] Weaver | Diffuse fields in open systems and the emergence of the Green's function(L)[END_REF]; [START_REF] Sanchez-Sesma | Retrieval of the Green's function from cross correlation: the canonical elastic problem[END_REF]; [START_REF] Wapenaar | Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF]; [START_REF] Roux | Extraction of coherent wavefronts from ocean ambient noise[END_REF]). However, these optimal conditions are almost never fulfilled in seismology as the noise sources (primarily oceanic, [START_REF] Gerstoft | A year of microseisms in southern California[END_REF]; [START_REF] Landès | Origin of deep ocean microseisms by using teleseismic body waves[END_REF]) are distributed unevenly. The consequence is that the Green function is only partially reconstructed, and that a certain amount of unpredicted fluctuations remain in the long-time averaged noise correlation. However, using dense network of seismic stations, the large number of correlation pairs provides enough redundancy from the data to counterbalance the relative inaccuracy of traveltime measurements extracted from noise correlation. As a matter of fact, seismic tomography based on noise correlation has been performed successfully at various scales ranging from kilometers to hundred of kilometers [START_REF] Shapiro | Highresolution surface-wave tomography from ambient seismic noise[END_REF]). Another application of seismic noise correlation is the monitoring of seismic velocity variations. Indeed, an advantage of using noise correlation to retrieve the Green function is that it is a continuous measurement. This allows to track temporal changes from a set of correlation functions calculated on successive time windows of noise recordings. For example, by looking at small phase changes in the signals extracted from noise correlation, we may detect velocity changes occurring in the medium sampled by the seismic waves. In practice, the coda portion of the noise correlation is used for measuring seismic velocity changes. The coda waves have been scattered throughout the medium, and have accumulated the delays due to a velocity change along their trajectories. This makes them more suitable to measuring velocity variations than the direct waves. What makes this use of noise correlation even more appealing is the work performed by [START_REF] Hadziioannou | Stability of Monitoring Weak Changes in Multiply Scattering Media with Ambient Noise Correlation: Laboratory Experiments[END_REF] through laboratory experiments that showed that passive monitoring was successful without completely reconstructing the GF from cross-correlations. This means that, contrarily to passive imaging, passive monitoring would be viable with a few stations and an uneven distribution of noise sources. Actually, two ingredients are necessary for passive monitoring. First, measuring accurate phase variations between similar signals requires a good signal-to-noise (SNR) ratio. In particular, the velocity variations associated with earthquakes are very subtle -of the order of 0.1% relative change -except for the non-linear response of unconsolidated materials at the surface [START_REF] Peng | Temporal changes of shallow seismic velocity around the Karadere-Düzce branch of the north Anatolian fault and strong ground motion[END_REF]; [START_REF] Rubinstein | Seismic velocity reductions caused by the 2003 Tokachi-Oki earthquake[END_REF]; [START_REF] Niu | Migration of seismic scatterers associated with the 1993 Parkfield aseismic transient event[END_REF]) or in highly damaged fault zones [START_REF] Li | Seismic velocity variations on the San Andreas fault caused by the 2004 M6 Parkfield Earthquake and their implications[END_REF]). Since the coda waves strongly contribute to the accuracy of phase variation measurements, the challenge is therefore to obtain a good SNR on this part of the correlation function that is classically of much lower amplitudes than the direct waves. Second, since the spatial distribution of seismic noise is uneven, the correlation function may evolve with time as the noise sources are modified (Froment et al. [START_REF] Larose | Fluctuations of correlations and Green's function reconstruction: role of scattering[END_REF]; [START_REF] Sabra | Surface wave tomography from microseisms in southern california[END_REF]). This need for long averages might be problematic though, since we would like to follow short term variations in the medium as well. In particular, Brenguier et al. (2008a) have measured the velocity changes occurring around an array near Parkfield, California, in the period spanning from 2003-2007. They observed drops in seismic velocity around the array which seem to coincide with the 2004 M w = 6.0 Parkfield earthquake, as well as with the 2003 M w = 6.5 San Simeon earthquake. However, since the noise was averaged over 30 days to retrieve the Green function, it is very difficult to say if the observed changes are coseismic or slightly delayed or advanced.

As for the time evolution of the noise source distribution, one may claim that the coda part of the noise correlation function will be less sensitive to the spatial distribution of the noise sources than the direct waves, as was experimentally shown at the laboratory scale [START_REF] Hadziioannou | Stability of Monitoring Weak Changes in Multiply Scattering Media with Ambient Noise Correlation: Laboratory Experiments[END_REF]).

In conclusion, the key issue for passive monitoring is the combination of a strong SNR on strongly-scattered waves extracted from noise correlation.

To reach such condition, different timescales play a role in the success of passive monitoring which will be discussed throughout the paper (see table 1): T SNR is the averaging time of the noise correlation function, T medium is the time evolution of the medium properties (the local velocity change after an earthquake, for example) and T noise is the time evolution of the spatial distribution of the noise sources.

In this paper, we revisit the 2003-2007 Parkfield data introducing a new processing algorithm with the objective to improve the temporal resolution of the velocity change measurements. The goal of the processing is to improve the SNR of our signals while keeping short temporal averages. The paper is structured as follows: in the next section (section 2), we discuss the data set and the signal processing methods previously applied on the noise data. Then, (section 3), we introduce an adaptive filter to improve the SNR of noise-extracted correlation functions, which increases the temporal resolution of velocity variation measurements (sections 4 and 5). We then investigate if the origin of the velocity variations can be attributed to changes in the ambient noise source distribution (section 6). Finally, in section 7 we introduce an expression for estimating the errors made on velocity variation measurements. Conclusion is that the velocity drop observed with waves in the band from 0.1 to 0.9 Hz at the Parkfield earthquake is coseismic, and that the distribution of the noise source directions has no influence on the velocity variations observed when using the coda of the signal.

Parkfield Data

The data used in this study is recorded at the Berkeley High Resolution Seismic Network (HRSN3 ) near Parkfield, California. This network consists of 13 borehole stations which have been recording continuously since July 2001. The borehole depths range from 60 m to 300 m, which helps reduce surface effects such as temperature variations and precipitation. The location of the stations and these two earthquakes is shown in Figure 1. For each pair of stations, we consider the vertical component of the continuous record between the 1st of August 2003 andDecember 31st, 2006. In this period, the M W = 6.5 San Simeon earthquake occurred on 22 December 2003, as well as the M W = 6.0 Parkfield earthquake of September 28, 2004. The noise signal recorded at each station is whitened between [0.1 0.9] Hz using a tapered boxcar window in the frequency domain, and subsequently 1-bit digitized. The daily cross-correlation is then computed for each of the 78 station pairs as well as the 13 autocorrelations, for lapse times between -120 and +120 seconds.

Initially, in an effort to improve the SNR, the daily cross-correlations are stacked with a sliding window of 30 days, which is shifted 1 day at a time. Later on, sliding window stacks of 1 and 5 days are used as well. An example of the signals obtained for a 30 day sliding window for one station pair (JCNB-SMNB) is shown in Figure 2. For each station pair, we define a reference correlation as the average of the correlations for the whole 1248 day measurement period. The temporal evolution of the seismic velocity can be measured by comparing the coda of the reference correlation to that of the shorter stacks. Temporal variations are measured using two different methods, as detailed in section 4. Seasonal variations are visible on the ballistic waves (Figure 2, center) of the cross-correlation as an asymmetric change in intensity. This suggests a seasonal change in noise source, as will be investigated in section 6. Note The gray horizontal dashed lines represent the dates of the San Simeon earthquake (22 December,2003) and the Parkfield earthquake (28 September, 2004). In the center a zoom on the direct waves, on either side on the coda part of the correlations. Note that the seasonal variation visible on the direct waves is no longer present in the coda. that this intensity variation is not visible in the coda part of the signal (Figure 2, left and right).

Method: Adaptive Filter

In an attempt to improve the SNR of the noise correlations, thereby increasing the temporal resolution, an adaptive filter (developed by [START_REF] Baig | Denoising seismic noise cross correlations[END_REF]) is applied to the correlations before retrieving the velocity variations.

The S-transform

The algorithm used is based on the S-transform, as developed by Stockwell [START_REF] Stockwell | Localization of the complex spectrum: The S transform[END_REF]). The analytical S-transform of a continuous function of time x(t), noted as S{x}(τ, f ), is defined as:

S{x}(τ, f ) = 1 √ 2π ∞ -∞ x(t)|f | exp - (τ -t) 2 f 2 2 -i2πf t dt, (1)
which transforms the signal x(t) into a time-frequency domain (τ, f ). The Gaussian envelope function (exp[-(τ -t) 2 f 2 /2]) localizes in time, and is translated along the signal. In the meantime, the oscillatory exponential kernel exp[-i2πf t] fixes the frequency being analyzed, and does not translate. This allows the S-transform to find the amplitude and phase independently. In other words, the S-transform not only estimates the local power spectrum, but also the local phase spectrum. For a discrete, real signal x(t) with N points, the S-transform will result in N 2 /2 complex-valued points (the positive frequencies). The S-transforms therefore needs to be evaluated at N 2 /2 points in time-frequency space. [START_REF] Stockwell | A basis for efficient representation of the S-transform[END_REF] developed the Discrete Orthonormal S-transform (DOST) to reduce the number of evaluations necessary. The N -point signal is projected onto a series of N complex basis functions x k [τ, f ]; k = 1 : N , with each function mapping to a localized patch of time-frequency space described by its frequency, f , and a time lag, τ . These basis functions form an orthonormal set. Like its parent S-transform, DOST allows both straightforward manipulation of signals from time to time-frequency while preserving the phase of the signal.

Both the analytical S-transform and the DOST are linear, which is beneficial in application to noise discrimination. Modelling a recorded time series, x(t) as signal, g(t), corrupted with noise, n(t), (ie. x(t) = g(t) + n(t)), the Stransform of the signal (S{x}(τ, f )) becomes:

S{x}(τ, f ) = S{g}(τ, f ) + S{n}(τ, f ) (2)
Provided the signal and noise occupy different portions of time-frequency space, the absence of cross-terms in equation 2 potentially allows for the noise cancellation. We shall exploit this property in section 3.2.

Construction of the Self-Filter

To construct a denoising filter based on the data itself we use the linearity property of the S-transform, as shown in equation 2. The challenge is to identify the signal g(t) when only the noisy timeseries x i (t) are measured. Suppose we have M records of the same timeseries, x i (t). When we average over all these x i (t), and assume a random noise, the noise part of the record should tend to 0 as the number of records approaches ∞, and only the signal remains. We can see that with equation (2) the S-transform of the timeseries becomes:

S{x}(τ, f ) = S{g}(τ, f ) + 1 M M k=1 S{n k }(τ, f ) (3)
where . represents a normalized sum over the signals 1

M M i=1 .
The second term in equation 3 disappears as M → ∞.

We can now construct a phase coherence filter F , inspired by the method used in [START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF], using the 'cleanest' possible signal x(t). The signal x(t) is decomposed into its DOST basis vectors x k . The normalized, complex values of these basis vectors are stacked:

F (τ, f ) = 1 N N k=1 S{x k }(τ, f ) |S{x k }(τ, f )| ν , (4)
The parameter ν can be used to control the strength of the filter. In this paper, ν = 0.5. In this case, the reference stack takes the role of the clean signal x(t) used in equation ( 3) and is used to construct the filter for each station pair. Each daily cross-correlation is then S-transformed, filtered according to:

∞ -∞ S{x i }(τ, f )F (τ, f )dτ = X i (f ), (5)
and finally transformed back to the time-domain.

A drawback of the analytical S-transform (equation 1) is that it is computationally expensive to calculate. In practice, we apply the DOST transform to our signals in order to gain on computing time.

Method: Doublets & Stretching

Two processing techniques have been proposed in the literature to estimate relative velocity changes dv/v in the coda of the correlations. The first one, the doublet technique, was developed by [START_REF] Poupinet | Monitoring Velocity Variations in the Crust Using Earthquake Doublets -An Application to the Calaveras Fault, California[END_REF]. It is also known as cross-spectral moving-window technique (CSMWT, [START_REF] Frechet | Application of the Cross-Spectral Moving-Window Technique (CSMWT) to the Seismic Monitoring of Forced Fluid Migration in a Rock Mass[END_REF]), and in practice computes the phase shift between records for consecutive, overlapping time windows ∆t. For a given window ∆t, the time shift δt is assumed to be constant and is estimated in the frequency domain by measuring the Fourier cross-spectrum phase. This estimator uses an accurate, unbiased Wiener filter technique [START_REF] Jenkins | Spectral Analysis and Its Applications[END_REF]) and produces an estimate whose confidence interval is controlled by the coherence values in the frequency band used for the analysis. The method can then measure arbitrary time-shifts between two records with enough similarity (or coherence). It is consequently possible to represent the time-shift δt versus time t in the coda. The actual velocity variation is simply the average slope of δt(t): dv/v = -δt/t. The key parameter in this analysis is the Fourier transform window length ∆t, which is the time-window over which the time-delay is measured. On one hand, large time windows include more data, thus averaging down the fluctuations of δt due to noise. This results in more accurate measurements. On the other hand, the approximation that the time-shift is constant within the given time-window is getting more erroneous when the window length is increased. The choice of ∆t is therefore a trade-off between time-shift accuracy, and the time resolution between two consecutive measurements. This processing found remarkable applications in seismology, including recent developments in ambient noise seismology (Brenguier et al. (2008b,a)).

Another idea [START_REF] Lobkis | Coda-Wave Interferometry in Finite Solids: Recovery of P-to-S Conversion Rates in an Elastodynamic Billiard[END_REF]; Sens-Schönfelder and Wegler ( 2006)) is to interpolate the coda at times t(1 -ε) with various stretching factors ε.

The actual relative velocity change is the stretching factor ε 0 = dv/v that maximizes the cross-correlation coefficient:

X(ε) = t 2 t 1 h [t(1 -ε)] h 0 [t]dt t 2 t 1 h 2 [t(1 -ε)] dt. t 2 t 1 h 2 0 [t]dt , (6)
where t 1 and t 2 the start and end time of the coda used, respectively. One drawback is that this latter processing assumes a linear behavior δt versus t, equivalent to a homogeneous relative velocity change dv/v = ε. This is sometimes not the case in complex media. But this technique has also a noticeable advantage: whole data is processed at once, which is found to result in a more stable, and thus more precise, estimation of dv/v. In particular, we showed from laboratory experiments that the stretching technique is more adapted to data with low signal-to-noise ratio [START_REF] Hadziioannou | Stability of Monitoring Weak Changes in Multiply Scattering Media with Ambient Noise Correlation: Laboratory Experiments[END_REF]). It was also found very useful for noisy seismic correlations on the moon (Sens-Schönfelder and Larose ( 2008)). Another interesting feature of the stretching technique is to provide the remnant coherence X(ε). This coherence indicates if noise sources are stable or changing over the period of observation of interest: X ∼ 1 means absolutely stable sources, X << 1 means that sources locations are changing. X also indicates the quality of the ε estimation (see section 7). In the case of Parkfield, we end up with a dv/v resolution of 10 -3 for measurements on 1 day stacks, and up to 2 • 10 -4 for 30 day stacks (section 7).

Results

The adaptive filter is tested with the two methods (doublets and stretching) for retrieving velocity variations detailed in section 4. Both analyses are done on cross-correlations bandpassed for [0.1 0.9]Hz.

Velocity Variations from Doublets

In the doublet analysis the time shift δt is measured for 100 windows of ∆t = 10s, which are shifted by 1 second. A slope is then fitted along the measurements for time windows ranging from -120s to -20s and +20s to +120s in the cross-correlation signal. The window from -20s to 20s is excluded, in order to minimize effects from the direct waves. Figure 4: dv/v using the doublet method, for a stack of 30 days. In blue: using the original data, in red: after applying the adaptive filter. The dashed line represents displacement measurements along the fault measured at the pomm GPS station (figure 1). The first vertical dashed line indicates the date of the San Simeon earthquake, the second that of the Parkfield earthquake.

The extent of the stack used is shown as a shaded gray area.

Figure 4 shows the relative velocity changes obtained with the doublet method on 30 day stacks (which slide by 1 day), before (blue) and after (red) application of the adaptive filter. We can see that while the dispersion of the measurements is reduced, as seen by the reduction of the errorbars, the overall values of velocity change measured are reduced. This follows from the way the doublet method retrieves the time lag in the frequency domain [START_REF] Poupinet | Monitoring Velocity Variations in the Crust Using Earthquake Doublets -An Application to the Calaveras Fault, California[END_REF]). The delay of each small window of the signal is measured by fitting a line to the phase-frequency behavior in the cross-spectrum of the signal and its reference. In Figure 3, the absolute value of each of the N DOST basis vectors is portrayed on its corresponding time-frequency space. In a way, it represents the absolute value of the DOST basis function for the reference signal of station pair JCNB-SMNB.

In this figure, it is visible that the adaptive filter applies different filtering coefficients for different frequencies, which affects this linear phase-frequency behavior. Since the doublet method measures the time delay by fitting a linear regression to the phase-frequency, this way of filtering interferes with the measurements. This lower amplitude of the dv/v is a disadvantage of the adaptive filter when it is used in combination with the doublet technique. However, in some cases it can be useful to use the filter to increase the SNR and thereby the temporal resolution. This way, it is possible to pinpoint the date at which a change occurred, even if the amplitude is wrong. For stacks shorter than 10 days, the SNR of the signals become too small to apply the doublet method, even after filtering. Figure 5: dv/v using the stretching method, for a stack of 30 days (A), 5 days (B) and 1 day (C). In blue: using the original data, in red: after applying the adaptive filter. The black dashed line represents displacement measurements along the fault measured at the pomm GPS station (figure 1). The first black dashed line indicates the date of the San Simeon earthquake, the second that of the Parkfield earthquake. The shaded gray areas around the dashed lines show the extent of the stack used.

Velocity Variations from Stretching

The stretching technique measures the dilation of the signal using as much of the coda as possible (see equation 6) -in this case, from ± 20s to ± 120 s. This renders the technique less sensitive to fluctuations in the correlations.

In Figure 5a we show the velocity changes measured with the stretching method using 30 day stacks, with and without application of the filter. Again, we remark a reduction of fluctuations after applying the filter. However, unlike in the case of the doublet method, the filter does not degrade the amplitude of the variations measured. In Figure 5b, we observe the same with a stack of 5 days. Note that the postseismic behavior is still consistent with the postseismic slip as measured by [START_REF] Freed | Afterslip (and only afterslip) following the 2004 Parkfield, California, earthquake[END_REF] from GPS measurements4 . Another improvement by the filter can be seen when re- Figure 6: dv/v using the stretching method, for a stack of 1 day, zoom around Parkfield earthquake (28 September, 2004). In blue: using the original data, in red: after applying the adaptive filter. The gray shading represents the day of the earthquake. Due to the many aftershocks this day interfering with the background noise, the dv/v measurement is unreliable.

The points for the day of the earthquake are therefore not plotted. The errorbars are described in section 7.

ducing the stack length to only 1 day (Figure 5c). Here, when using the unfiltered signal, the velocity drops associated with the San Simeon and Parkfield earthquakes can be observed, but they become much clearer in the case of the filtered signal. This is especially apparent when we look closely at the velocity behavior around the date of the Parkfield earthquake (Figure 6). The errorbars shown in this figure are described in section 7.

The measurements from the filtered signals show that the velocity one day after the earthquake has dropped to a lower value, and recovers slowly.

Note that the ambient noise on the day of the earthquake is altered by numerous aftershocks and the mainshock itself. After elimination of these parts of the noise, not much data remains to work with on that day, thus producing an unreliable correlation signal. There also appears to be a velocity increase one day before the event, however, we would be very hesitant to claim a precursor based on this single data point. From these observations (Figure 6), we can safely conclude that the velocity drop we measure has its maximum at the time of the Parkfield earthquake [START_REF] Rubinstein | Depth constraints on nonlinear strong ground motion from the 2004 Parkfield earthquake[END_REF]).

Stretching, especially in combination with the adaptive filter to remove as much noise from the signals as possible, allows to greatly improve the temporal resolution of the measurements.

Noise Source Location from Beamforming

After the previous discussion on the improvement of the temporal resolution of the dv/v measurements, we study if the observed fluctuations in the velocity can be explained by the change in noise source distribution. We observe on the example of cross-correlations of couple JCNB-SMNB (Figure 2, center) that the direct waves are sensitive to seasonal variation, probably due to variations in the noise source distribution. According to the theory [START_REF] Gouédard | Convergence of the two-point correlation function toward the green's function in the context of a seismic-prospecting data set[END_REF], [START_REF] Froment | Estimation of the effect of non-isotropic distributed energy on the apparent arrival time in correlations[END_REF], [START_REF] Paul | Empirical synthesis of time-asymmetrical Green functions from the correlation of coda waves[END_REF]), we expect any strongly directive and varying part of the noise field to influence the Green function reconstruction negatively and thus alter the dv/v measurement.

In order to establish if any such strong moving sources exist around the Parkfield array, we investigate the structure and evolution of the ambient noise field. To do so, we perform plane wave beamforming (as described in Roux ( 2009)) on 1-day segments of ambient noise in the [0.1 0.2]Hz frequency band which contains the secondary microseism peak. Figure 7 shows an example of the angular-speed distribution of the background noise on this frequency interval. We observe that the ambient noise field at Parkfield is very directive and dominated by a component propagating to ∼ 50 • clockwise from North with an apparent velocity around 2500 m/s. This means that it consists mostly of surface waves originating from the direction of the Pacific Ocean, in agreement with results from similar analyses by [START_REF] Roux | Passive seismic imaging with directive ambient noise: application to surface waves and the san andreas fault in parkfield, ca[END_REF].

We track the azimuthal variation of the most prominent maximum in Figure 8a. In order to quantify the temporal change of the incident direction of all noise sources, not only the strongest, we define a correlation of the noise incident direction. This is evaluated by looking at the daily coherence of the beamformer output with respect to the output averaged over the three years of data (Figure 8b). Do the azimuthal variations of the dominant source, or variations of the noise field coherence, reflect on the dv/v measurements? In Figure 8, we have plotted the azimuthal variation of the strong source at [0.1 0.2]Hz, as well as the noise field coherence (a and b) along with the dv/v measurements and the correlation function coherence (c and d).

From Figure 8a, some seasonal variations in the azimuth of the beamformer output is visible, reminiscent of the variations seen on the cross-correlations in Figure 2. However, one can see that the change in noise direction does not influence the dv/v measurements directly. On the other hand, the change in noise field structure, expressed as a decoherence in Figure 8b, does relate to a change in correlation waveform (Figure 8d). As seen in section 4, this decrease in waveform coherence as the source distribution changes, is expected. In some cases, such as in the beginning of 2005, a strong change in noise sources has even led to more fluctuations in the dv/v measurements. This effect of decreased waveform coherence on the measurements is included in the expected rms fluctuations, as can be seen in section 7.

In conclusion, a change in noise wavefield incidence does temporarily result in a noisier dv/v measurement. Since such a change is also visible as a decoherence of the waveform these periods of time can be detected and treated accordingly.

Error calculation

The precision with which a dilation coefficient ǫ is retrieved depends on the strength of correlation between two signals. [START_REF] Weaver | On the precision of noise-correlation interferometry[END_REF] derived an expression to predict the fluctuations in dilation coefficient due to e.g. changes in the noise source, in the case where the velocity in the medium has not changed. The goal is to evaluate the precision with which wave speed changes can be evaluated. To do this we consider the case in which two waveforms differ only by noise so that the actual relative dilation, without noise, is zero. We then look for the apparent (non-zero in general) value of ǫ at which the corresponding X in equation ( 6) achieves its maximum. Equation 7gives the root mean square of this apparent, and erroneous, relative dilation.

(ǫ) = √ 1 -X 2 2X 6 π 2 T ω 2 c (t 3 2 -t 3 1 ) , ( 7
)
where ǫ is the dilation, X the coherence of the waveform after dilation, between starting time t 1 and end time t 2 . T is the half period bandwidth, with ω c the center frequency of the signal. In this study, values are typically ω c = 0.5, T = 0.4, with start time t 1 = 20 seconds and end time t 2 = 120 seconds. This leads to an rms(ǫ) evaluated by: rms

(ǫ) = 2.6 × 10 -3 √ 1 -X 2 2X , (8)
as detailed in [START_REF] Weaver | On the precision of noise-correlation interferometry[END_REF]. Since we average the measurement over 91 station pairs, this rms value is reduced by √ 91, which leads to a rms of the order of 8 × 10 -5 for a wavefield coherence of X = 0.8. If the dilation measured exceeds this rms value, we can conclude that it is indeed due to a velocity change in the medium. We find that in the case of the Parkfield data, the velocity variations measured are larger than the predicted fluctuation value of ǫ for a correlation coefficient of X. This means that they are due to a local velocity change and not loss of coherence in the cross-correlation function.

Discussion

By improving the temporal resolution of our dv/v measurements using an adaptive filter and the stretching technique, we have shown that the velocity drop observed at Parkfield is coseismic. The increase in temporal resolution obtained in this paper allows us to track velocity changes immediately after the event. Similar, albeit stronger, coseismic velocity drops have been observed and are generally associated with shallow damage in either the crust [START_REF] Rubinstein | Evidence for Widespread Nonlinear Strong Ground Motion in the MW 6.9 Loma Prieta Earthquake[END_REF], [START_REF] Sawazaki | Time-Lapse Changes of Seismic Velocity in the Shallow Ground Caused by Strong Ground Motion Shock of the 2000 Western-Tottori Earthquake, Japan, as Revealed from Coda Deconvolution Analysis[END_REF], [START_REF] Peng | Temporal changes of shallow seismic velocity around the Karadere-Düzce branch of the north Anatolian fault and strong ground motion[END_REF], [START_REF] Niu | Migration of seismic scatterers associated with the 1993 Parkfield aseismic transient event[END_REF]) or the fault zone [START_REF] Li | Postseismic Fault Healing on the Rupture Zone of the 1999 M 7.1 Hector Mine, California[END_REF]). These shallow effects are modeled by a strong drop in shear modulus (and thus velocity drop) very shortly after the mainshock, followed by a postseismic logarithmic healing curve. The peak ground acceleration measured in the region surrounding the array is in the order of 0.65 -1.3g [START_REF] Shakal | Preliminary analysis of strong-motion recordings from the 28 September 2004 Parkfield, California earthquake[END_REF]), so superficial damage is a possiblity. It could be argued that the velocity drop we observe is a result of this shallow damage. The lower amplitude of velocity drop -on the order of 0.1% instead of 30% -could be the result of averaging a localized effect over a larger volume, seeing as we measure the change with multiply scattered coda waves. However, if the velocity drop we observe were due only to superficial damage, we would expect the drop amplitude to become higher as we measure in a shorter time window just following the event (1 day as opposed to 30 days). As seen in Figure 5, this is not the case. In fact, the postseismic relaxation behaviour closely follows along-fault displacement deduced from GPS measurements [START_REF] Freed | Afterslip (and only afterslip) following the 2004 Parkfield, California, earthquake[END_REF], [START_REF] Johanson | Coseismic and Postseismic Slip of the 2004 Parkfield Earthquake from Space-Geodetic Data[END_REF]). This suggests the possibility that not shallow, but deeper, deformationrelated effects are involved here. Apart from the velocity drop at the Parkfield event, some small velocity fluctuations remain. As shown in section 7, these fluctuations are larger than the error expected on dv/v, and should have a physical origin. Since our GF is not perfectly reconstructed, these fluctuations could be due to a change of the sources. In section 6, we find no direct relation between the dominant noise direction and the dv/v measurements (Figure 8). However, the amount of decorrelation of the cross-correlations does reflect the stability of the ambient noise wavefield.

Conclusion

By using an adaptive filter to rid the data of as much noise as possible, and by applying the stretching method to retrieve the velocity variation measurements, we are able to perform monitoring on cross-correlations of only 1 day of ambient noise. Thanks to this, we can show that the 0.1% velocity change observed at the Parkfield event is well above the measurement precision and precisely coseismic. Furthermore, we show that slow seasonal variations in the orientation of the main noise sources do not affect the velocity change measurements directly. Sudden changes in the noise wavefield do relate to large fluctuations measured in dv/v.

Figure 1 :Figure 2 :

 12 Figure 1: Location of the HRSN stations (white and black triangles) near Parkfield, California. In the inset figure, the locations of the 2003 San Simeon and 2004 Parkfield earthquakes (stars). The white and black circle indicates the location of the GPS station pomm.

FilterFigure 3 :

 3 Figure 3: The absolute value of the coefficients of DOST basis vectors for the reference signal for JCNB-SMNB with each number coefficient plotted at its appropriate patch of time-frequency space. The colorscale ranges from dark blue -representing zero amplitude -to dark red for the highest amplitudes.

Figure 7 :

 7 Figure 7: Beamforming output representing the angle-versus-velocity ambiguity surface at [0.1 0.2]Hz. Azimuth in degrees clockwise from North.

Figure 8 :

 8 Figure 8: (a) The azimuth of the beamformer output maximum (in degrees clockwise from North) and (b) the beamformer output coherence, both at [0.1 0.2]Hz. (c) The velocity variations as well as (d) the correlation waveform coherence for a 5-day sliding stack. The gray vertical dashed lines represent the dates of the San Simeon earthquake and that of the Parkfield earthquake, respectively.

  rms

Table 1 :

 1 Overview of the different timescales used in this paper. For example, seasonal changes of ocean noise directivity have been observed, whose effect on the direct arrivals of the noise correlation function is comparable or of greater amplitude than an actual change of the medium velocity. Going back to the SNR issue, improving SNR in the seismic noise correlation function means averaging over longer noise time series of typically weeks or months

	designation typical duration definition
	T	0.1 -1 s	Dominant period of the correlated signal
	t	10 -120 s	Lapse time of the correlated signal
	T SNR	1 -30 days	Averaging time of cross-correlation; stack length.
	T noise	days -months	Evolution timescale of the noise sources
	T medium	months -years	Evolution timescale of crustal properties
	(2010)).

Institut des Sciences de la Terre, Université Joseph Fourier and CNRS, BP53, 38041 Grenoble Cedex 9, France

Now at: ESG Solutions, 20 Hyperion Court, Kingston, ON, K7K 7K2, Canada

More information concerning the Parkfield HRSN can be found at http://seismo.berkeley.edu/bdsn/hrsn.overview.html

More information concerning the USGS GPS network at Parkfield can be found at http://earthquake.usgs.gov/monitoring/edm/parkfield/continuous.php

Acknowledgments

This work was funded by the ANR JC08 313906 SISDIF grant, and by the ERC advanced grant WHISPER. The authors thank F. Brenguier, P. Gouédard, N. Schaeffer, R. Weaver and N. Shapiro for fruitful discussions and useful remarks. All the seismic data used in this study came from the Parkfield HRSN and were collected by the Berkeley Seismological Laboratory (BSL).


Diffuse transport in clay media: µm to nm scale characterization of pore space and mineral spatial organization S. Gaboreau 1 , J.C. Robinet 2 , Ch. Tournassat 1 , S. Savoye 3

1. BRGM (French Geological Survey), 3, Avenue Claude Guillemin, 45060 Orleans Cedex, France 2. ANDRA, 1-7 Rue Jean Monnet, 92298 Châtenay-Malabry, France 3 CEA, L3MR, CEA/Saclay, 91191 Gif sur Yvette, France

In the framework of radioactive waste repository, clayrock formations are foreseen as barrier materials due to their diffusion properties. In clay materials, the dominant transport mode is diffusive and depends mainly on various parameters such as the mobility of the species in water, the accessible porosity, the pore space geometry and the retardation as a result of reactions such as sorption or ion exchange (Tournassat and Appelo, 2011). In this way, the European CATCLAY project (EURATOM FP7), in the context with research on transport in porous materials, was proposed to describe the cation migration processes in natural clayrocks. The project is structured along 3 RTD workpackages, combining modeling and experimental studies from a simpler, analogous system (monophasic compacted clay system) to clayrocks (Callovo-Oxfordian argillites, Opalinus Clay and Boom Clay).

Part of this experimental studies focuses on small scale structure (µm -nm) property of rocks in order to determine how the spatial distribution of mineral and pores at small scales can influence diffusion driven transport of sorbing cations. The present study focuses on compacted illite properties (simpler analogous system) in hopes to extent this study to the natural clayrock formation. Illite was chosen by the way that is the main constituent of clayrock. Compacted illite material represents thus an analogy with the clay matrix constituting clay-rocks. Our approach is mainly based on imaging the small scale structural organization of compacted illite material and analyzing the obtained images in order to extract information on pore space and mineral spatial distribution.

Techniques for imaging the texture of illite material like water saturated, in compacted state, were first developed. The first step was to improve classic resin impregnation method in order to preserve the texture without losing the clay confinement and modifying the pore space geometry. This has been done by taking into account the molecule size of the monomer, the low viscosity, the dipole moment (adapted for the clayrock with swelling clay content) and the controlled time polymerization. MMA monomer proved to be the most suitable resin in our study. The small scale structure of impregnated sample was then imaged in 2D using Transmission Electron Microscopy (TEM) and in 3D using Focused Ion Beam coupled to Scanning Electron Microscopy (FIB/SEM). For TEM observations, a set of ultra-thin serial sections (50 -100 nm) were cut using a microtome. A set of 2D images were then acquired using a resolution ranged between 100 nm and 10 Å. TEM images clearly show us the multi-scale organization of clay materials (Figure 1 and 2); we observe the 10 Å spacing sheets constituting the illite particles, nanometer size illite clay particles more or less aggregated and the surrounding pores having a size ranging from few hundred nanometers to nanometer.

FIB/SEM analysis is currently in progress. From FIB/SEM, a set of serial images can be acquired using the "slice and view" method (Keller et al., 2011). Then, 2D FIB/SEM images need to be aligned to reconstruct a 3D volume. Image resolution is limited to 10-20 nm. Both methodologies (FIBtomography and TEM techniques) are thus complementary method for the up-scaling characterization of the structural organization of compacted clayey materials. TEM images analysis allow to scale down the resolution size since only a part of the pore space could thus be imaged with FIB/SEM method (Keller et al., 2011).

Viewing and performing a qualitative description of images constitute a major result and can help us to better understand how the transfer pathways and retention sites are organized in the porous media.

Thanks to image analysis method, pores and minerals can be thresholded from grey level TEM and

  Figure 1.

  Figure 2.


Alain Dutech
email: alain.dutech@loria.fr

Self-organizing developmental reinforcement learning

This paper presents a developmental reinforcement learning framework aimed at exploring rich, complex and large sensorimotor spaces. The core of this architecture is made of a function approximator based on a Dynamic Self-Organizing Map (DSOM). The life-long online learning property of the DSOM allows us to take a developmental approach to learning a robotic task: the perception and motor skills of the robot can grow in richness and complexity during learning. This architecture is tested on a robotic task that looks simple but is still challenging for reinforcement learning.

Overview

The framework of reinforcement learning [START_REF] Sutton | Reinforcement Learning[END_REF]] is particularly attractive because it enables an artificial agent to learn an action plan to solve a problem that may be unfamiliar or uncertain while only using a scalar signal to distinguish between good and bad situations. Learning is not supervised and does not require an expert or an oracle. Moreover, Markov Decision Processes (MDP) [START_REF] Puterman | Markov Decision Processes: discrete stochastic dynamic programming[END_REF] provide a formal background for analyzing theoretical properties of many algorithms, ensuring, for example, the existence of an optimal solution to the learning problem and convergence to this optimal solution.

Although reinforcement learning led to some successful applications, the practical use of reinforcement learning for realistic problems is not easy. Difficulties usually arise from the size of the problems to be solved that can only be described through a very large number of states, making them computationally intractable for reinforcement learning algorithms searching an exact solution. Algorithms providing approximate solutions do exist (approximate value iteration, LSPI, direct policy search through gradient ascent, etc., see Chap. 3,4 of [START_REF] Sigaud | Markov Decision Processes in Artificial Intelligence[END_REF], [START_REF] Lagoudakis | Least-squares methods in reinforcement learning for control[END_REF]), but they are far from providing answers to all problems encountered in practice.

We are particularly interested in the use of reinforcement learning in the context of autonomous robotics. Several difficulties reduce the applicability of the learning algorithms. The sensor data are continuous in value, often noisy, expensive to obtain (in time and energy) and very (too) rich in information (e.g. in the case of a video stream). Our main concern is for the agent to be able to explore efficiently these rich and complex environments. We propose a developmental approach where the agent starts with crude actions and perceptions that gradually get more complex. Thus, starting from a smaller and simpler sensorimotor space, the agent should be able to learn simple tasks and, building over acquired behaviors, learn more complex tasks as its sensorimotor world gets richer. This developmental learning is made possible by using a self-organizing adaptive map architecture for function approximation in a reinforcement learning framework.

This paper first quickly presents the reinforcement learning framework used and original architecture for a developmental approach (Sections 2 and 3). Then, the experimental setting is detailed and results are given (Sections 4 and 5). Section 6 discusses the our work, pointing out some limitations and offering perspectives for future work.

Principles and Motivations

Our approach of developmental reinforcement learning is largely based on the Q-learning algorithm [START_REF] Watkins | Learning from delayed rewards[END_REF] and therefore relies on estimating the optimal value function Q * in the sensorimotor space S × A. The principle of the algorithm is to use each training sample (s t , a t , s t+1 , r t ) provided at time t through an interaction between the agent and its environment (the agent perceives s t , chooses action a t that changes the environment now perceived as s t+1 , and gets a reward r t ) to update the estimate of the value function Q * for the pair (s t , a t ) as follows:

Q * t+1 (s t , a t ) ←-Q * t (s t , a t ) + ∆ α r + γ max a ∈A Q * t (s t+1 , a ) -Q * t (s t , a t ) , (1)
where α is a learning coefficient, which may depend on s, a and t.

A major limitation to the practical use of reinforcement learning methods is related to the complexity of algorithms. Theoretically, to guarantee the convergence of reinforcement learning algorithms to an optimal solution one requires that each (s, a) pair is visited an infinite number of times. Even with only a very large number of visits, this constraint is one of the main limitations to the practical use of reinforcement learning techniques, and particularly in the context of robotics, and the only way to learn good behavior would be to explore intelligently and efficiently the environment.

The sensorimotor space of a robot is rich, continuous and complicated, furthermore it can be time and energy consuming to get a new training sample, and also risky (like bumping in a wall). As a result, an efficient and complete exploration of the unknown environment is very unlikely. The main idea of our developmental reinforcement learning framework is to ease the exploration by starting with a very small and limited sensorimotor space and gradually increasing it when the robot learned performances improve. Here, this support to the agent mainly takes the form of a gradual increase of perceptions richness and potential actions of the agent as its performance in the learning task is improving.

From this perspective, several problems arise. First, the learning architecture chosen must allow to progressively increase the sensorimotor space of the agent and the difficulty of the task. Second, we want the agent to rely on what he has already learned to learn in a richer environment. This last issue is particularly delicate, it is similar to the problem of knowledge transfer, a research field still wide open [START_REF] Caruana | Learning to Learn, chapter Multitask Learning[END_REF]. Within the framework of reinforcement learning, a central element of learning is the estimation of the value function of the sensorimotor space. The problems we have mentioned are related to estimating a function: how to take advantage of the current estimate of a function to extend it to a modified input space (when sensorimotor space is enriched).

Architecture for developmental learning

The core of our architecture is an adaptive function approximation to learn the Q-value function in a continuous sensorimotor space. The principle of our architecture is depicted in Fig. 1. Perceptions of the robot (we give more details on these perceptions in Section 4), denoted s, feed a Dynamic Self-Organizing Map (DSOM) [START_REF] Rougier | Dynamic Self-Organising Map[END_REF]. The activity of this self-organizing map is the input of a perceptron with one layer of neurons that has as many outputs as possible actions. Thus, the output neuron associated with the action a learns the value function Q(s, a) in a supervised way using the Q-learning update equation (1) a linear regression pattern, by modifying its weights ω using: ω ← ω -L ∆. The DSOM map uses unsupervised learning to adapt to the perceptual inputs received.

Although we have not used a more conventional architecture like LSTD or LSQ [START_REF] Bradtke | Linear least-squares algorithms for temporal difference learning[END_REF][START_REF] Lagoudakis | Least-squares methods in reinforcement learning for control[END_REF], our system is actually quite close. One can consider that the use of a DSOM map projects the perceptions on basis function Φ(s) that evolve over time, and that the linear regression estimates the value function as a linear combination w T Φ(s) of these basis functions. The search for optimal coefficients of this linear combination is made by a stochastic gradient descent and the basis functions being adaptable, our algorithm does not guarantee convergence. The neuronal learning architecture is well adapted to the developmental approach that we want to experience. The action space is discrete, its richness is closely tied to its cardinality. Increasing the richness of the action space can be dealt with by adding neurons to the output layer. These new neurons take advantage of the already adapted DSOM map. The learning of the coefficients of the linear combination of a new output neuron can be accelerated by initialising these coefficients using the coefficient of neurons coding "semantically" close actions. In the example that we present below (see Section 4), a new neuron connected with action turn-right-slowly will be initialized with the coefficients related to the neuron for the action turn-right.

For the continuous state space, the richness comes in part from to the dimension of this space. Increasing the size of the input space is not as simple as for the output space. One solution is to increase the size of all vectors of input weight without changing their projection on the previous perceptual space. Another solution that we would like to test is to provide the architecture with input weight of the maximum dimension of the input space and, while the current perceptual vector has dimension less than this maximum, to artificially increase its size by cloning a few of these values to the extra dimensions. The latter approach seems better adapted to transfer what is learned in small size to larger tasks.

Experimental platform

Our experiments are performed on a KheperaIII robot of K-Team company1 . Although this robot embarks a processor, we chose to deport processing and computation on an external computer, communication between this computer and the robot is done by wifi.

The robot moves using two drive wheels which allow it to turn on the spot. The robot can use discrete actions, such as moving forward or backward a certain distance (on average, as the actions of the robots are noisy) or turn right or left for a given angle. We use discrete actions because we want to explore a framework based on Q-Learning and estimating a value function for continuous actions is still largely an open problem. The robot's action space is a discrete space and its richness is characterized by the number of available actions.

The main sensor we use is a wireless camera placed above the robot which captures a 320x200 bitmap color image, the image is then processed to provide sensory information in a lower dimensional space. Through a logical sensor we call "retina", the robot is provided with perceptions made of a vector of dimension dim per where each coefficient of the vector, between 0 and 1, is the ratio of a particular hue computed on a vertical strip of the image. The number and positions of the vertical stripes are customizable and can evolve during the experiments. On the example of the right side of Fig. 2, the robot is facing a blue diamond on a red background, the retina is configured to show the ratio of blue vertical stripes -named B1, B2 and B3 -positioned at the abscissa 80, 160 and 240 of the image. Perceptions of the robot in this case are [0.42, 0.67, 0.33]. The state space of the robot is a continuous space with a richness measured by the dimension of this space (that is to say the number of vertical stripes). Given our experimental platform, the tasks we want to teach our robot are navigation tasks in mazes with visual cues. Our idea is to set visual cues in the environment to guide the robot, e.g. arrows that tell in which direction to head in order to reach the exit. If the robot needs a fine perception to distinguish and recognize the arrows, it does not need this degree of perception in order to be "attracted" by the highly contrasted symbols. In addition, the robot can learn "satisfactory" behaviors with only crude actions. Then, these behaviors can be refined and precised as actions and perceptions become richer. We believe that this kind of scenario, which can be broken down into sub-tasks of varying complexity, is well suited to the experimentation of the concepts of developmental reinforcement learning.

Experimental results

The learning architecture proposed is tested in a very simple testbed. This experiment, which only involves an increase in the richness of the action space of the robot, is only a preliminary to other experiments. Increase in perceptions space and in the complexity of the task are scheduled.

In this experiment, the robot is positioned in a closed environment. The surrounding walls are red and a blue diamond (the target) is set on one of the walls. The robot must choose between 5 discrete actions: turn-right and turn-left (about 34 degree turn), stop, turn-right-slowly and turn-left-slowly (about 20 degree). These action are noisy by nature because of the imprecision in the robot command. The perceptual space is made of 3 vertical stripes that return the ratio of blue. These stripes, called B1, B2 and B3 are respectively set to the extreme left, at the middle and at the extreme right of the image given by the camera of the robot, their relative angle to the front of the robot are -16, 0 and 16 degree. The robot has to learn to face the target and it receives a reward if it stops while the target is in front of it. If 0.7B1 + B2 + 0.7B3 > 1.4 and the robot stops, then it gets a reward of +1, otherwise the reward is zero. After getting a non-zero reward, of after 12 moves without success, the robot is repositioned randomly.

Learning samples come from real movements and perceptions on the robotic platform but our learning algorithm needs too many iterations and the sample set is not enough for that. As explained later, 10, 000 to 15, 000 iterations are needed for the algorithm to learn. Thus, we use, and more importantly reuse, the same set of 2, 000 samples generated on the robot. Each iteration of the algorithm is done with one sample randomly chosen from the sample set, with respect to the learning scenario (e.g., a sample must use one of the actions currently allowed to the robot).

Three different learning scenarii are tested. In the first one, the robot can only use 3 actions turn-left, turn-left and stop. In the second scenario, the robot can choose among the 5 different actions. In the "developmental " scenario, the robot starts learning with 3 actions during N B iterations then it can also use the 2 other actions, with some optional adaptation of learning parameters.

Results presented were obtained using the following parameters: 64 neurons for the DSOM map, with a "small world" neighborhood and at least on neighbor The agent stops when facing the target (toward the (1,1,1) vertex) and uses big turns when the target is far from the center (weak values for B1 or B3).

for every neuron (nb neur = 64, nb link = 1) ; D = 0.5 and η = 1 are the learning parameters of the DSOM map; α = 0.1 and γ = 0.9 are parameters for the reinforcement learning framework; the exploration policy used is an epsilongreedy policy, with π = 0.25. The learning parameter L of the linear regressor varies with the experiments (usually between 0.001 and 0.5). This parameter seems to have a more important influence on the quality of the results and on the speed with which performances are reached. These parameters have been set by hand, they give good results, as suggested by the learned policy depicted on Fig. 3.

For every experiment, the quality of learning is evaluated every 1, 000 iterations by testing the greedy policy on a fixed set of starting positions. The set of positions is defined so as to discriminate as much as possible the different policies. These positions are given by the following angles: 0, 2, 4,6,8,10,12,14,16,18,20,22,24,25,30,35,40,45,50,90,180 and their symmetric. The mean of the reward received from these positions is a measure of the quality of a policy, and thus of the learning. The best hand-made policy using 3 actions has a performance of 0.138 and, with 5 actions, one can reach a performance of 0.39. Fig. 4 shows the performance reached by the algorithm as a function of the number of iterations when using 3 or 5 actions, with L = 0.05 and L = 0.1. Each graph is made of 10 learning trajectories. The points are the performances measured and the line are sliding mean of these trajectories. For this set of parameters, as for most of the other settings tested, one can notice a large variability in the performance of the algorithm along a trajectory. The variability can be limited by using very low L value.

Number of Iterations

We want to point out that the learning speed is largely increased with the DSOM based architecture. In a previous work using a multi-layer perceptron and eligibility traces [START_REF] Sarzyniec | Apprentissage par renforcement développemental en robotique autonome[END_REF], 100, 000 iterations were needed before reaching good performances, whereas here only 5, 000 to 15, 000 are required.

The main focus of the experiments is the exploration of the developmental scenario. Fig. 5 compares the performances of the algorithm when using directly 5 actions with performances obtained with the developmental scenario (3 then 5 actions). The direct learning is quicker but the developmental learning can rather quickly (around 10, 000 iterations) reach better performances. To reach these performances, the learning parameter L of the linear regressor must be adapted during the developmental scenario: starting at a value of 0.1, it is reduced to 0.01 when adding the two new actions. If this parameter is kept at 0.1, the performance are reduced (red curve in figure 5, right).

Fig. 5:right shows that, in a developmental scenario, the number if iterations before adding the last 2 actions does not influence the reached performance. For this set of parameters ( L = 0.1 and π = 0.25), as for many others, performances tend to decrease after 20, 000 iterations. This fact, combined with the large variability within a learning trajectory, may suggest that the task to learn is not that simple. This may be related to the difficulties of using non-linear approximators in a reinforcement learning framework [START_REF] Bertsekas | Neuro-dynamic programming[END_REF][START_REF] Coulom | Reinforcement learning using Neural Networks, with Applications to Motor Control[END_REF].

Discussion

Results presented here are still preliminary and are part of a work in progress. Many parameters influence the learning algorithm and these parameters should be more systematically explored. Some of these parameters seem to balance each other, like, for example, the number of neurons in the DSOM map and its learning parameters: with a large number of neurons, the DSOM map does not need to be very adaptive and D can be low. The intrinsic properties of DSOM map, like the capacity to adapt continuously to the input while not being oversensible to the distribution of the input, seem particularly interesting when the perceptual space of the agent will change with time. We want to test this, by increasing the perceptions space (add more stripes to the retina) but also by moving these stripes or making them more noisy.

Learning performances are very unstable. In only 100 iterations, the performance of the learned policy can change a lot. Even when the DSOM map is stabilized, which happens very quickly in our experiments as the perceptual space is stable and "small", performances are still unstable, even for low learning parameters of the linear regressor. We would like to study more systematically how the various parameters could be tuned to decrease this variability while preserving both the level of performance and the learning speed, but these criteria are qualitatively assessed right now. We lack an automated procedure for comparing the overall performance (speed, stability, level) of different parameter settings. This kind of automated overall performance evaluation could also be used to guide and control the increase of the richness of the sensorimotor space of the agent. Thus, using original developmental path, it could be possible to link the increase in performance and the "maturation" process of the agent, as done in [START_REF] Baranes | Maturationally-constrained competence-based intrinsically motivated learning[END_REF] for example.

Conclusion

In this paper, we have presented the principles of a developmental reinforcement learning architecture in order to ease the exploration of large and complicated sensorimotor spaces. At the heart of this architecture lies a dynamic self-organizing map that participates in learning an approximation of the value function of the sensorimotor space. Thus, our learning algorithm is a kind of developmental neuro Q-Learning algorithm.

Experiments conducted so far on a simple robotic task are quite promising but far from complete, a lot of testing is still needed. When increasing the action space, the learning architecture learns more quickly than previous architectures [START_REF] Sarzyniec | Apprentissage par renforcement développemental en robotique autonome[END_REF][START_REF] Dutech | Dynamic reservoir for developmental reinforcement learning[END_REF]. Furthermore, even on this simple task, the developmental scenario tested brings better results than a direct approach. This work opens a lot of perspective and future work, in particular a more systematic study of the parameters, the testing of evolving perceptual spaces and more complex tasks.

Figure 1 :

 1 Figure 1: Learning architecture. Robot perceptions s are the input of a dynamic self-organizing map. A linear output layer with one neuron for each actiona learns Q(s, a).

Figure 2 :

 2 Figure2: Left: the robot in its environment. Right: the perception given by the retina of the robot. For every vertical stripe, positioned at (80, 160, 240), the ratio of blue is computed, here [0, 42; 0, 67; 0, 33].

Figure 3 :

 3 Figure 3: Representation of a policy learned with 5 actions in the sensorimotor space (the coordinate B1, B2 and B3 are the dimension of the perceptual space).The agent stops when facing the target (toward the (1,1,1) vertex) and uses big turns when the target is far from the center (weak values for B1 or B3).

Figure 4 :

 4 Figure 4: Learning with 3 and 5 actions, points depict the performances every 100 iterations, for 10 different trajectories. Each trajectory is smoothed and represented by a line.

Figure 5 :

 5 Figure 5: Developmental learning. Left: The direct (5 actions) and developmental learning (3 actions, then 2 more after N B = 5000 iterations, L goes from 0.1 to 0.01) are compared. The curve show the mean and variance for 10 trajectories. Right: Developmental learning with various delay before using the full action set ( L = 0.1, N B = {5000, 10000, 15000}.

See http://www.k-team.com.


Emmanuel Promayon

Vincent Luboz

Grégory Chagnon

Thierry Alonso

Denis Favier
email: dfavier@grenoble-inp.fr

Christine Barthod
email: christine.barthod@univ-savoie.fr

Yohan Payan
email: ypayan@imag.fr

COMPARISON OF LASTIC (LIGHT ASPIRATION DEVICE FOR IN VIVO SOFT TISSUE CHARACTERIZATION) WITH CLASSIC TENSILE TESTS

LASTIC is a device estimating in vivo soft tissue elasticity. It uses negative pressure to deform the tissue surface and captures several deformation stages to trace the behavioral curve. Using Finite Element inverse analysis and a Neo Hookean constitutive law, the tissue's Young modulus is evaluated. This paper compares LASTIC capabilities with standard tensile tests on four samples with elastic properties ranging from 10 kPa to 100 kPa. Although LASTIC overestimates Young modulus by an average of 24 %, it allows a first estimation of the elastic modulus of different materials.

1.

INTRODUCTION

Modeling living soft tissues and their complex specificities (non-linear, inhomogeneous, anisotropic, patient-specific behaviors) is a challenging task, especially during a surgical intervention. The Light Aspiration device for in vivo Soft TIssue Characterization (LASTIC) was, in its first version, used to evaluate the constitutive behavior of skin, tongue and brain tissues [1]. It is based on the pipette aspiration principle: after putting the instrument in contact with the tissues surface, negative pressures are applied and displacement responses are measured. Consequently, the device had to be (1) fully sterilizable and (2) sufficiently fast to provide an interactive estimation of the tissues elastic parameters. These two constraints lead to the new version of LASTIC [2]. It is validated here by comparing it, for in-vitro materials, with usual tensile test machines. This validation is a prerequisite before any intra-operative use of this new version.

METHODS

Elastic samples

Experiments were conducted on four different samples with elasticity moduli ranging from few kilopascals to several hundred kilopascals, to simulate the range of human soft tissues' elasticity. RTV#1: a RTV-EC00 silicone (artificina.com), made from a mix of 50 % of base and 50 % of catalyst (also called curing agent). It has a linear behavior up to 25 % of engineering strain. RTV#2: the same RTV-EC00 silicone, but made from a mix of 40 % of base and 60 % of catalyst, to create a softer silicone than RTV#1. It has a linear behavior up to 25 % of engineering strain. Ecoflex: an Ecoflex® 00-30 silicone (Smooth-on.com) constituted by two bases mixed in equal proportion. It has a linear behavior up to 15 % of engineering strain. Candle gel: a gel (glorex.com) that is fluid at temperatures over 95 °C. It has a linear behavior up to 10 % of engineering strain. Two samples were created for each material in order to perform two different tests: one with LASTIC and one with the tensile tests, see subsections 2.2 and 2.3 for details on the mold shapes. The same material mix was poured in the two molds so that the elasticity measurements were done on the same material for each testing device.

LASTIC

The current version of LASTIC is designed to be used in vivo and consequently can undergo sterilization. LASTIC is built in a compact metallic cylinder of 33 mm in height and 34 mm in diameter divided in two compartments, see Figure 1. The lower one is airtight, open at the bottom by a circular aperture and closed at the top by a glass window. The upper compartment holds a 9.5 mm × 9.5 mm 2 megapixels camera (model VS6750, STMicroelectronics) with a resolution of 1600 × 1200 pixels. The negative pressure ∆P is created in the lower compartment by a programmable syringe pump (Aladdin AL-1000, World Precision Instruments), while a USB digital manometer (model 8215, AZ Instrument Corp) measures the pressure. As a result, the soft-tissue is aspired inside the chamber, where the deformation is imaged by the digital camera, by means of a 45° inclined mirror. Increasing levels of negative pressure are applied during the measurement, creating larger aspiration height. The deformation, measured as the aspiration height, is segmented on each image of the measurement. It is converted into mm using a basic camera calibration determining the pixel size. A Finite Element Analysis of the aspiration experiment using a Neo-Hookean constitutive law [3] is used in order to build a library of displacement heights [4]. The library of displacements depends on Young modulus and the applied negative pressure. From this library, a least-square minimization method is used to find, in less than 1 s, the Young modulus in the library that best matches the measurements.

The samples used for LASTIC tests are molded as cylindrical shapes of approximately 20 mm of height, 60 mm of diameter.

3.

RESULTS

The results for LASTIC and the tensile devices are given in details in [4] for each material. Figure 1 (right panel) summarizes the Young modulus estimated by each device. For the candle gel, Young moduli of 10.5 kPa and 14 kPa were found for the tensile devices and LASTIC respectively, resulting in 33 % difference. Similarly, for the other materials, differences of 28 %, 22 % and 11 % were found for the RTV#2, Ecoflex, and RTV#1 respectively. On average, LASTIC thus overestimates the Young modulus by 23.8 % ± 9.5 % SD compared to the tensile devices.

4.

DISCUSSIONS

The average discrepancy of 23.8 % in the Young modulus E estimation is relatively high. However, this overestimation seems to be linearly dependent of the value of E and this could be taken into account to improve soft tissue characterization of living tissues with LASTIC. This linear variation could be caused by our digital manometer. It has indeed a precision of ± 6 mbar in measuring the negative pressure created with the pump. Although rather small, this still represents 3 % of the maximum applied negative pressure for the RTV#1, 5 % for the Ecoflex, 7 % for the RTV#2, and 15 % for the candle gel. Since these variations are higher for more elastic materials, our assumption that the manometer is responsible for part of the measurement deviation could be corroborated. A more accurate manometer could help to reduce this error. Nevertheless, for a precise estimation of the Young modulus, we should also improve LASTIC on two other sources of error: the camera calibration and the choice of the constitutive law for the inverse problem solving as explained in [4]. Finally, it is important to keep in mind that, while the tensile tests give the closest approximation available of the elastic modulus of the different materials, they are not gold standards. For example, the estimation of the stress/strain curve of each material with a Neo Hookean constitutive equation is not perfect as the measurements are not exactly reproducible.

5.

CONCLUSIONS

Overall, despites an average overestimation of 24 % compared to standard tensile devices, LASTIC gives interesting estimations of the Young modulus of elastic materials by aspirating an area of their surface and recording the corresponding displacement. This device could already be used on patients or even in surgical set up because of the overestimation consistency and because of LASTIC's small size, small cost and capability to be sterilized. Future improvements, presented in the discussions section, should provide accurate in vivo soft tissues measurements.

Figure 1 -

 1 Figure 1 -LASTIC cross section: the upper part contains the camera while the lower part is the aspiration chamber with the mirror (left). Young moduli found by LASTIC and the tensile devices for each material (right). 2.3 Tensile devices Two extension devices were used for tensile measurements to compare with LASTIC tests. Both apply a given displacement to the material, fixed between two tensile jaws, and record the force needed to reach it. The SYMME Laboratory uses an Instron® machine (Instron, Norwood, MA, USA). The material sample had a rectangular shape of 60 mm x 15 mm x 7 mm. The 3S-R laboratory uses an Eplexor 500N testing device (Gabo Qualiter Testanlagen GmbH, Ahlden/Aller, Germany), with a 25 N sensor. The material had a rectangular shape of 25 mm x 7 mm x 2 mm. To evaluate the elasticity modulus from the measurements of both devices, the tensile test stress/strain values are studied by fitting a Neo-Hookean constitutive equation [4].

  Light Aspiration Device. Medical Image Analysis, 2009, 13:673-678. 2. Schiavone, P., Promayon, E. and Payan, Y. LASTIC: A Light Aspiration Device for in vivo Soft TIssue Characterization, Biomedical Simulation: 5th International Symposium, ISBMS, 2010, 5958:1-10. 3. Treloar, L. R. G. The elasticity of a network of long chains molecules I and II.Trans. Faraday Soc., 1943, 39, 236-246. 4. Luboz V., Promayon E., Chagnon G. et al. Validation of a Light Aspiration device for in vivo Soft TIssue Characterization (LASTIC). Soft Tissue Biomechanical Modeling for Computer Assisted Surgery, Springer-Verlag, 2012.

	Material	Young modulus in kPa dif in %
		Tensile	LASTIC
	Candle gel	10.5	14	+33
	RTV#2	25	32	+28
	Ecoflex	55	67.5	+22
	RTV#1	90	100	+11
	6.

REFERENCES

1. Schiavone, P., Chassat, F., Boudou, T., Promayon, E., Valvidia, F. and Payan, Y. In Vivo Measurement of Human Brain Elasticity Using a