Coda Wave Interferometry for Estimating Nonlinear Behavior in Seismic Velocity

See allHide authors and affiliations

Science  22 Mar 2002:
Vol. 295, Issue 5563, pp. 2253-2255
DOI: 10.1126/science.1070015


In coda wave interferometry, one records multiply scattered waves at a limited number of receivers to infer changes in the medium over time. With this technique, we have determined the nonlinear dependence of the seismic velocity in granite on temperature and the associated acoustic emissions. This technique can be used in warning mode, to detect the presence of temporal changes in the medium, or in diagnostic mode, where the temporal change in the medium is quantified.

In many applications, such as nondestructive testing or monitoring of volcanoes or radioactive waste disposal sites, one is primarily interested in detecting temporal changes in the structure of the medium. Temporal changes in Earth's structure that accompany earthquakes have been observed on the basis of the attenuation of coda waves (1), on the arrival times of the directly arriving waves (2), on velocity changes inferred from later arriving waves (3) [see also (4)], and on changes in seismic anisotropy (5). Here, we introduce coda wave interferometry whereby multiply scattered waves are used to detect temporal changes in a medium by using the scattering medium as an interferometer. For quasi-random perturbations of the positions of point scatterers, or for a change in the source location or the wave velocity, estimates of this perturbation can be derived from multiply scattered waves by a cross correlation in the time domain.

In the numerical example (Fig. 1), the wave field for a medium consisting of isotropic point scatterers is computed with the use of a deterministic variant (6, 7) of Foldy's method (8). Given the mean free path (l = 20.1 m) and the wave velocity (v = 1500 m/s), one can infer that aftert = 5.4 × 10−2 s the waves are on average scattered more than three times. The later part of the signal is called the coda. Suppose that one repeats this multiple scattering experiment after the scatterer locations are perturbed. The perturbation in the scatterer location is 1/30 of the dominant wavelength and is uncorrelated between scatterers (9).

Figure 1

Location of 100 scatterers before and after the perturbation (filled circles and open circles, respectively) with the source (asterisk) and receiver location (triangle). For the sake of clarity, the scatterer displacement is exaggerated by a factor 40. The scatterers are placed in an area of 40 m by 80 m. The waveforms recorded before and after the perturbation at the receiver are shown on the right with a solid and dashed line, respectively.

In this example, the scatterers' locations are perturbed. In general, a perturbation can involve other changes in the medium or a change in source location. We refer to the waveform before the perturbation as the unperturbed signal and to the waveform after the perturbation as the perturbed signal. For early times (t < 0.04 s), the waves in Fig. 1 have not scattered often, rendering the path lengths of these waves insensitive to the small perturbations of the scatterers (small compared with the dominant wavelength λ = 2.5 m), which causes the unperturbed and perturbed signals to be similar. However, the multiply scattered waves are increasingly sensitive with time to the perturbations of the scatterer locations because the waves bounce more often among scatterers as time increases. The correlation between the unperturbed and perturbed signals, therefore, decreases with increasing time.

The perturbation in the medium can be retrieved from the cross correlation of the coda waves recorded before and after the perturbation. The unperturbed wave fieldu unp(t) can be written as a Feynman path summation (10) over all possible paths P Embedded Image(1)where a path is defined as a sequence of scatterers that is encountered, tP is the travel time along path P, AP is the corresponding amplitude, and S(t) is the source wavelet. When the perturbation of the scatterer locations (or source location) is much smaller than the mean free path, the effect of this perturbation on the geometrical spreading and the scattering strength can be ignored, and the dominant effect on the waveform arises from the change in the travel time τP of the wave that travels along each pathEmbedded Image(2)The time-windowed correlation coefficient is computed fromEmbedded Image Embedded Image(3)where the time window is centered at time t with duration 2T andt s is the time shift used in the cross correlation. When Eqs. 1 and 2 are inserted, double sums ΣPP over all paths appear. In these double sums, the cross terms with different paths (PP′) are incoherent and average out to zero when the mean of the source signal vanishes. This means that in this approximationEmbedded Image(4)where ΣP(t,T) denotes a sum over the paths with arrival times within the time window of the cross correlation, and the autocorrelation of the source signal is defined asC(t) ≡ ∫−∞ S(t′+t)S(t′) dt ′.

For time shifts τ much smaller than the dominant period, a second-order Taylor expansion givesC(τ) = C(0)(1 − 1/2ω̄2 τ2), whereω̄2 is the mean-squared frequency of the multiply scattered waves that arrive in the time window. Using this givesEmbedded Image(5)where <… >(t,T) denotes the average for the wave paths with arrivals in the time interval (tT,t+T).

The time-shifted cross correlationR (t,T)(t s) has a maximum whenEmbedded Image(6)where 〈τ 〉(t,T)is the mean travel time perturbation of the arrivals in the time window. The value of the cross correlation at its maximum is given byEmbedded Image(7)where στ 2 is the variance of the travel time perturbations for waves arriving within the time window. Therefore, the mean and the variance of the travel time perturbation of the waves arriving in the time window can be extracted from the data recorded with a repeatable source and one or more receivers.

Different types of perturbations leave a different imprint on the time-shifted correlation coefficient. When the scatterer locations are perturbed independently with root mean square displacement δ, the mean travel time perturbation vanishes (<τ >(t,T) = 0), and the variance is given by (7)Embedded Image(8)where l * is the transport mean free path (11). In deriving Eq. 8 we use that the number of scatterers encountered is on average given byn = vt/l, where t is the time that the wave has spent in the scattering medium. Using Eqs. 7and 8, the root mean square perturbation of the scatterer location follows from the maximum of the time-windowed correlation coefficientEmbedded Image(9)A different type of perturbation is a constant change δv in the velocity for fixed locations of the scatterers. The mean travel time perturbation is given by 〈τ 〉(t,T) = − (δv/v)t, and when the time window is small (Tt), στ = 0. The velocity change follows from the time of the maximum of the time-shifted cross correlation functionEmbedded Image(10)When the perturbation consists of a displacement of the source location over a distance δ for a fixed medium, only the wave path to the first scatterer is perturbed. In that case, the mean travel time perturbation vanishes 〈τ 〉(t,T) = 0, and for an isotropic source the variance is given by 〈στ 2〉 = (δ/v)2. The source displacement, then, follows fromEmbedded Image(11)These different perturbations can be distinguished on the basis of the time-shifted cross correlation. When the positions of the scatterers are perturbed, the mean travel time perturbation vanishes and the maximum of the cross correlation decreases linearly with increasing time, whereas for the perturbation of the source position the maximum value of this function is independent of time. A change in the velocity is detectable by a shift in the position of the maximum ofR (t,T)(t s), which increases linearly with time.

The root mean square displacement of the scatterers inferred from the numerical example (Fig. 1) is shown in Fig. 2 as a function of the center timet of the time window. The inferred change δ in the scatterer location does not depend on the center time of the window used for the cross correlation. This provides a consistency check of the method.

Figure 2

The value of δ obtained from the time-windowed cross correlation of the waveforms in Fig. 1 and of 20 other receivers as a function of the center time t of the time window (solid line), plus or minus one standard deviation (dotted lines) for T = 2 × 10−2 s. The true root mean square displacement value δtrue = 8 × 10−2 m is shown by the horizontal solid line.

The extreme sensitivity of the coda waves to changes in the medium is used here in a laboratory experiment to infer the temperature dependence of the seismic velocity in Elberton granite. In many experiments, the change in the seismic velocity in rock samples is measured for a temperature change of about 100°C (12, 13). In our experiment, a cylindrical sample of granite with a height of 110 mm and a diameter of 55 mm was heated from 20° to 90°C with a heating coil inside the sample and then was cooled down to room temperature. The heating and cooling phase each took about 8 hours. Two piezo-electric transducers were used to excite and record elastic waves in the sample with a dominant frequency of about 100 kHz. The waveforms were recorded after each ±5°C change in temperature. In order to reduce the influence of ambient noise, the waveforms were averaged over 10 repeated measurements. A third transducer was used to monitor the acoustic emissions in the sample.

The difference in the early part of the waveforms recorded at temperatures of 45° and 50°C (Fig. 3) are small. This change in temperature does not affect the first arrival, which means that the travel time of the first arrival cannot be used to infer any possible small change in velocity due to a 5°C temperature difference. The late time window (Fig. 3, bottom inset) shows a clear time shift of the waveforms.

Figure 3

Waveforms recorded in the granite sample for temperatures of 45° and 50°C, in blue and in red, respectively. The insets show details of the waveforms around the first arrival (top inset) and in the late coda (bottom inset).

For each change of ±5°C in temperature, the change in the velocity is inferred from Eq. 10 using 20 different time windows of the coda waves, each with a duration of 0.1 ms. The mean and variance of the velocity change (Fig. 4) is inferred from the estimates of the velocity change in the different time windows. The relative velocity change is of the order of 0.1% for a temperature change of ±5°C with an error of about 0.02%.

Figure 4

The absolute value of the relative velocity change for a 5°C increase and 5°C decrease, red and blue symbols, respectively, as a function of the highest temperature during the change. The histograms shows the count of acoustic emissions.

During the heating phase, the velocity change is constant for temperatures below 75°C. Above that temperature, the velocity change increases during heating (Fig. 4). The acoustic emissions correlate with the increased value of the velocity change at 75°C (14). During the cooling phase, the velocity change is constant and there are no acoustic emissions. When the sample is heated again to a temperature of 90°C, the velocity change does not increase dramatically around 75°C and there are no acoustic emissions (15). In order to test whether the transducer coupling and the presence of the heating coil played a role, we repeated the experiment with an aluminum sample. In that case the velocity change is constant both during heating and cooling.

The acoustic emissions and the change in the velocity gradient occur only in a pristine sample during heating [the Kaiser effect (14)] and are due to the irreversible formation of fractures by differential thermal expansion (16) of the minerals in the sample. This indicates that the velocity change is due to two different mechanisms. The first is a reversible change in velocity due to the change in bulk elastic constants with temperature. The second mechanism is associated with irreversible changes in the sample that generate acoustic emissions. The damage done to the sample leads to a greater change in the seismic velocity with increasing temperature.

These measurements could be carried out because of the extreme sensitivity of coda wave interferometry to changes in the medium. This technique makes it possible to infer the nonlinear dependence of the velocity on temperature that is associated with irreversible damage to the granite sample.


View Abstract

Navigate This Article