An experimental study of phase angle ﬂuctuation in seismic waves in random heterogeneous media: time-series analysis based on multivariate AR model

SUMMARY Effects of small-scale heterogeneities on seismic waveform ﬂuctuations were studied by physical model experiments. Using a laser Doppler vibrometer, we recorded elastic waves propagating through a granite block at 180 observation points that were arranged as an equally spaced circular array. A disc-shaped PZT source was attached on the other side surface of the circular array for realizing equivalent positions with respect to both source radiation pattern and travel distances of waves. Waveform pairs were selected out from the 180 waveforms, and cross spectra of time-windowed partial waveforms were calculated by applying the multivariate AR model. By comparing waveforms of two observation points, the cross-spectral amplitudes and phases are obtained with respect to the lapse time by moving the time window, or to the spatial distance by changing the pairs of observation points. We obtain distributions of cross-spectral phase values for frequency and the lapse time of waveforms. The distributions indicate phase ﬂuctuation of waves in random media with respect to frequency and lapse time. Heterogeneity of the rock sample is expressed as a 1-D exponential autocorrelation functions (ACF); (cid:2) 2 exp( −| r | / a ), where r is the distance, and a and (cid:2) are the correlation length (0.22 mm) and the strength of heterogeneity (8.5 per cent), respectively. The distributions are plotted against ka ; the product of wavenumber and correlation length. For small ka , the distributions of phase are close to the Gaussian distributions with small variances, but the variances quickly become large above ka ≈ 0.2–0.3. Then the distributions become uniform between − π and π . This suggests that the incoherent scattered waves become dominant above a critical ka value (or a critical frequency for a medium), and phase information in later portions of waveforms will be lost. This may be important for extracting reﬂection, refraction or converted waves that are assigned as signals from geologic discontinuities because those signals may be strongly distorted by scattered waves produced from the small-scale heterogeneities of earth’s media.


I N T RO D U C T I O N
In subsurface imaging techniques, calculated waveforms are mostly based on large-scale heterogeneities (truncated model structures). However, real structures contain heterogeneities with scale-lengths smaller than the reconstructed model structure. Observed waves are, therefore, not exactly the same as those expected from the model structure because seismic waves fluctuate by scattering and diffraction from these small-scale heterogeneities (Aki & Richards 1980;Sato & Fehler 1998). For example, there are traveltime fluctuations in waves caused by the small-scale heterogeneities (Mukerji et al. 1995;Snieder & Aldridge 1995;Shapiro et al. 1996;Samuelides 1998;Sivaji et al. 2001a;Spetzler et al. 2002;Spetzler & Snieder 2004).
We then consider that subsurface structures are divided into two different scales of heterogeneities: one is the large-scale heterogeneity that can be determined by conventional seismic explorations and the other is the small-scale heterogeneity which is a smallscale variation from an estimated structure. The small-scale heterogeneities are considered random because their scale sizes are smaller than resolutions of conventional seismic methods. Only their statistical properties can be known. The media containing small-scale heterogeneities are, hereafter, referred to as random media.
When considering wave propagation in random media, there are two approaches: one comes from optical ray theory and the other comes from wave propagations in equivalent homogeneous media (Mukerji et al. 1995). Most of earlier works were based on the first approach where heterogeneity scale lengths are large compared to seismic wavelengths (Snieder & Aldridge 1995). Numerical simulations are often applied in the first approach. Recently, works based on the second approach have become common by employing laboratory physical model experiments (Sivaji et al. 2001a). Numerical methods have difficulties for the second approach because the number of grid points is too large for ordinary computers to implement calculations for small-scale random heterogeneities. Physical model experiments, on the other hand, are suited for the second approach because real materials such as rocks contain various scales of random heterogeneity.
In random media, traveltime fluctuations are caused by phase angle perturbations in seismic waves (Sato & Fehler 1998;Shapiro et al. 1996;Samuelides 1998). Waveform fluctuations can be recognized as variations of wave parameters. Even for equal sourcereceiver distances, amplitudes and phases show variations in random media (Sato & Fehler 1998;Sivaji et al. 2002).
Wave fluctuations can be clearly seen in the frequency domain as amplitude and phase variations are functions of frequency. Both fluctuations are described by distributions of amplitude spectra or phase spectra for waveforms observed at observation points equally spaced from a source. For example, if we assume a uniform distribution for phase fluctuation, fluctuations of phase in the waveform appear with an equal probability for phase angles between −π and π . In this case, the waves would be completely incoherent.
There can be many kinds of distributions that produce fluctuations in seismic waves (e.g. the Gaussian distributions of different variances). If we consider waveform pairs, fluctuations can be illustrated by distributions of the phase and amplitude cross spectra between waves. It is therefore necessary to investigate the probability distribution functions of phase and amplitude fluctuations in seismic waves, and to relate the distributions to the intensity and scale-length of random heterogeneities.
If we can reveal a relationship between wave fluctuations and parameters of small-scale random heterogeneities, we can apply this relationship to evaluating misfits between observed waves and the waves calculated from truncated model structures. Evaluating the misfits will be important for discussing the accuracies of model structures.
In random media, observed waveforms differ from each others even when observation points have an equal path-length in terms of an equivalent homogeneous medium and have an equal radiation pattern with respect to a source . Here we try to find a relationship between seismic wave phase fluctuations and random heterogeneity scale, by observing waveforms at several different points having the same source radiation pattern and the same source-receiver distance. A fine-grained granitic rock was used as a model for the random medium and seismic waves were detected by a laser Doppler vibrometer (Nishizawa et al. 1997).

Outline of laboratory experiments on seismic wave propagation in random media
To study relationships between seismic wave fluctuations and random heterogeneities, we have three approaches: theoretical studies, numerical simulations and laboratory experiments. Theoretical studies always give clear-cut results but often assume the ideal case that the whole medium can be described as the one containing a single type of heterogeneity. Numerical simulations can provide versa-tility applicable to most realistic cases where media contain different types of heterogeneity. However, it is still hard to perform numerical simulations for 3-D cases, because such simulations need huge computer resources. Laboratory experiments, on the other hand, can easily model 3-D cases. Most of experiments so far performed have given only qualitative results because the waveform measurements were not accurate enough to be used for quantitative analyses.
Recently, a laser Doppler vibrometer (LDV) was used for accurate measurements of elastic vibrations. LDV can be applied to laboratory model experiments for measuring elastic waves in the ultrasonic frequency range (Nishizawa et al. 1997(Nishizawa et al. , 1998Scales & van Wijk 1999). LDV measures the velocity of vibration of a solid surface within a very small area (less than 50 μm in diameter). No corrections are necessary for frequency responses up to 1 MHz (Nishizawa et al. 1998). LDV is also free from disturbances of wave fields because of its non-contact method for detecting surface vibrations. Others have applied LDV or similar methods to study seismic wave propagation in media containing random heterogeneities. (Sivaji et al. 2001a(Sivaji et al. ,b, 2002Fukushima et al. 2003;Scales & Malcolm 2003). Using LDV as a waveform detector, laboratory experiments are very promising for studying the effects of random heterogeneities on seismic wave fluctuations. We can use the waveforms in model experiments for quantitative analysis of seismic wave fluctuations in random media.

Random medium
When we use a model material for laboratory physical model experiments, we need to describe the random heterogeneity of the material. We selected Westerly granite as a model of a random medium. A technique for quantifying random heterogeneity in rocks was shown by Spetzler et al. (2002) and Sivaji et al. (2002). They obtained surface images of the rock samples, and converted the images to the mineral distributions by changing the original image colours into three-valued colours: white, grey and black that correspond to feldspar, quartz and biotite, respectively. Traverse lines of the mineral distribution give 1-D velocity distributions, and then they provide an average velocity fluctuation in the rock sample.
According to their analysis, the 1-D velocity fluctuation in Westerly granite can be expressed as the exponential autocorrelation function: where r, and a are the distance, the intensity of fluctuation and the correlation distance, respectively. and a are measured as 8.5 per cent and 0.22 mm, respectively. Random heterogeneity of Westerly granite and other related parameters are shown in Table 1.  Figure 1. An illustration of physical model experiment viewing a section of the sample block. Elastic waves are generated from a piezoelectric transducer (PZT) attached at one of the major surfaces of the sample and received at the opposite surface. Waveforms are measured over a circular array of 20 mm diameter. Observation points are located at every 2 • with its centre coinciding with the symmetric axis of the disc-shaped PZT. Expected ray paths for direct P wave and reflected P wave (PP P ) are shown.

Waveform measurements
We use the same experimental data used in Sivaji et al. (2002). The details of the waveform measurements are already described in that paper. Here we describe only a few important points. Fig. 1 illustrates source and receiver geometry. A disc-type piezoelectric transducer was attached at the centre of the square surface (300 × 300 mm) of a 90 mm thick rectangular rock prism. The transducer has a 5-mm diameter and 2-MHz resonance along the thickness direction, radiating P and S waves by producing the normal force at the interface (Tang et al. 1994). P and S wave radiations are axisymmetric with the unique axis normal to the surface of transducer. Elastic waves propagating through the rock sample were detected by LDV at the opposite surface of the transducer. The laser beam irradiates observation points on a 10-mm radius circle centred on the normal of the disc-type transducer. The circumferential interval of observation points was two degrees, which corresponds to a nearest-neighbour distance of about 0.3 mm.
Waveforms were sampled by a fast A/D converter with a sampling rate 20 ns. Original waveform data was 12-bit resolution with a sample length of 4096 points. Data acquisition at each observation point was repeated 2000 times. All waveform data were summed and averaged to remove high-frequency electric noise. Figs 2(a)-(d) show the waves at 180 observation points for different source signals that are one-or two-cycle sine waves of different frequencies. The plotted amplitudes are in arbitrary scales; the maximum peaks of waveforms range from 0.5 to 1 mm s −1 .
An optical unit of LDV was set on a tri-axial stage which was driven by stepping motors. The optical unit moves parallel to the sample surface along the two mutually orthogonal axes: horizontal and vertical directions. It also moves along the direction orthogonal to the sample surface for beam focusing. Locations of the optical unit were measured by linear encoders with an accuracy of 1 × 10 −6 m. We used a reflection sheet covered by small glass lenses to improve efficiency of beam reflection. The highest intensity of the laser beam was obtained when the beam was focused on the centre of a glass lens. The beam was not accurately located on the expected positions on the circle, because the position of the glass lens generally deviates from the expected position. The average diameter of the glass lenses is 50 × 10 −6 m, and the deviation of the laser beam from an expected position is roughly as same as the diameter of the lens. The deviation is a main source of positioning error of the laser beam and results in errors of phase measurements, depending on wavelengths of the elastic wave. Since Westerly granite has a P-wave velocity 4.78 km s −1 , the wavelength of the 2-MHz P wave is about 2.5 mm. Phase-detection error caused by positioning error is similar to ± (2πδr/λ), where ±δr is the positioning error, and λ is the wavelength. If δr is similar to the diameter of the glass lenses, the error amounts to ±0.04π radian. This value is negligible because we found larger phase difference between observed waveforms as will be shown later. Errors of the present measurements are more clearly seen by analysing waveforms in steel that is considered as a completely homogeneous material.
Figs 2(a)-(d) indicate waveform fluctuations when seismic waves are propagating through random media. The waves of 0.25-and 0.5-MHz sources show only small fluctuations and the reflected P wave (PP P ) is very clear. However, in the 1-and 2-MHz source waves, PP P is masked by strongly scattered waves dominated by almost the same frequencies as the source waves. PP P is produced by a large-scale structure of which dominant wavelength is similar to the thickness of the sample.
In all waveforms, there are incoherent waves due to the small scale heterogeneity in the rock sample. Those waves appear between the signal waves: P, S and PP P . Incoherent waves become stronger with increasing source-signal frequency.

1-D AR model
Before describing the multivariate AR model, we start with a 1-D AR model and explain the concept of Akaike Information Criterion, AIC. We express the 1-D AR (autoregressive) model as where M is the order of AR model, and v n is a white noise sequence that follows the normal distribution with mean 0 and a variance σ 2 : N(0, σ 2 ). The data (x 1 , . . . , x N ) is one set of realizations from the model given by eq.
So far, we have assumed that the order of the AR model, M, is given. The best estimate of the AR order can be determined by searching the minimum AIC with respect to M: The criterion is derived from entropy maximizing principle which gives the best fit between the estimated distribution and the real distribution, or minimizes the divergence between the two distributions on the basis of the Kullback-Leibler information (Akaike & Kitagawa 1999).

AR model with dimension
We use -dimensional (multivariate) AR model of the order M described below: where x n is a vector of dimensional time series at the time n defined by: (Kitagawa & Gersch 1996;Akaike & Kitagawa 1999), and v n is the multivariate Gaussian white noise sequence with dimension having zero mean and the covariance matrix W . We use E (·) to denote the expected value of the term shown in the parentheses. The expected values associated with v n are the followings: where the prime indicates the transpose, O is a vector having elements zero, and W is an × symmetric matrix; σ 2 jk = σ 2 k j and is positive-definite.
We obtain the following cross spectrum matrix P(f ): where p jk (f ) is a complex number when j = k. Normally, define the Fourier transform of the AR operator by: where I is the × unit matrix. Then we have the spectrum matrix as: where A * denotes the complex conjugate of A. We use α jk for the absolute value of p jk (f ).
where, R{ p jk ( f )} and I{ p jk ( f )} are the real and imaginary parts of the spectrum p jk (f ). The phase of p jk (f ) is given by φ ij (f ) is calculated between −π and π.

Numerical calculation
We picked waveform pairs from the 180 waveforms observed at different points in the circular array. This corresponds to the 2-D multivariate AR model: = 2. From those waveform pairs, the data within 6-μs time window are used for the time-series analysis. The time window is moved forward with a step of 1 μs. We have 180 partial waveform pairs at each lapse time from the onset of the P-wave to the time 10 μs after PP P arrival. When picking waveform pairs, we change the distance between observation points by skipping array points. Distances between two observation points are described by numbers of points counted from the reference point along the clockwise direction: the numbers 1, 2, . . . indicate the pairs in which one is reference point and the other is the next neighbour, the second neighbour, and so forth, respectively. Hereafter those numbers are referred to as the station distance. The waveforms in the time window contain 300 data points with the sampling rate 20 ns. The sampling rate corresponds to the Nyquist frequency of 25 MHz. Since LDV shows a flat frequency response up to 2 MHz, the spectra are numerically valid up to this frequency.
We used a software package TIMSAC (Akaike et al. 1979) provided from the Institute of Statistical Mathematics. The multivariate AR coefficient matrix was calculated on the basis of AIC (Sakamoto et al. 1986;Kitagawa & Gersch 1996;Kitagawa 1999). AIC was used to objectively determine the order of AR model, M, in eq. (7). We set the maximum value of M as 17; M ≤ 17. Actually, the minimum AIC appeared between 9 and 17. In spectrum analysis, the order of M controls the number of peaks in spectrum. The present analysis is limited in the low-frequency region, so that the difference in the numbers of order M does not affect the results. We obtained 180 sets of parameters corresponding to the waveform pairs for each time window, and calculated α jk and φ jk as functions of frequency, spacing of observation points, and lapse time.

Waveform in Westerly granite
Waveform gathers in Fig. 2 show the direct P, the direct S, and the reflected, PP P , phases. Hereafter, we refer to those phases as signal waves. Incoherent waves are apparent between those signal waves. They also appear in the signal-wave parts and are regarded as the noise that masks signal waves. Incoherent waves become stronger as the frequency of source signal increases. The waves excited by the 2-MHz source frequency (Fig. 2d) contain strong incoherent waves that mask the S and the PP P phases. Incoherent waves are the scattered waves due to small-scale heterogeneity in the rock sample. The strength of scattered waves depends on the frequency of source signal and the lapse time. We analyse waveforms in the frequency domain by calculating cross spectra of waveform pairs. If the two waves have a strong correlation, the cross spectrum between the waves will indicate this. To study the correlation between waves, we first investigate how the amplitude and phase of the cross spectrum change as functions of the lapse time in waveform and the distance between two observation points.

Amplitude spectrum
By changing the station distance, we obtain other sets of data that may suggest effects of station distance on waveform fluctuation. Since the dimension of the multivariate AR model is two, the nondiagonal elements of the spectral matrix P(f ), p 12 (f ) and p 21 (f ), show relationships between the waveform pairs. Amplitude (α 12 (f )) and phase (φ 12 (f )) of the cross spectra p 12 (f ) are plotted for different lapse time and for different station distance.
Figs 3(a) and (b) show amplitude and phase of the cross spectrum at the time windows 21 μs (P) and 58 μs (PP P ) for the 0.25-, 0.5-, 1-and 2-MHz source signals (from the top to the bottom), respectively. The waveform pairs were taken for the station distance 7, which corresponds to the actual distance 2.4 mm. The cross spectra are shown in a frequency range from 0.125 to 3 MHz. All the spectra of 180 waveform pairs are plotted. Since the amplitude spectrum is plotted in an arbitrary absolute logarithmic scale, only relative intensities of spectra are significant.
Each α 12 (f ) has a peak that corresponds to the frequency of the source signal except for the case of the 2-MHz source signal, where α 12 (f ) (Figs 3(a.1-A) and (b.1-A) ) shows a peak at about 1.5 MHz, which is lower than the source-signal frequency. This is probably due to the strong scattering attenuation of the frequency components higher than 2-MHz. We hypothesize that this strong attenuation is due to the presence of thin microcracks. The existence of thin microcracks is suggested by a remarkable increase of P-wave velocity in Westerly granite under confining pressure (Nur & Simmons 1969). Since we have no LDV system equipped in a pressure vessel, we have no opportunity to confirm the conjecture about scattering from cracks by measuring waveforms under confining pressure. Even if we have this inconsistency of spectral peak in the measurements for a 2-MHz source, this is irrelevant to the present analyses because microcracks affect the wave fluctuation only at high frequencies.
The characteristic correlation length a and the intensity , which are obtained from the surface image, may be modified when we consider the effects from microcracks. However, the present measurements do not have an accuracy at frequencies higher than 2 MHz because of the positioning error of the laser beam (shown later in the measurements using a steel block). We therefore consider that the present experimental data are reliable up to 1.5 MHz.

Phase spectrum
φ 12 (f ) was calculated for the range (−π, π). The values shift abruptly from one side to the other side in the high frequency parts (Fig. 3) because of the cyclic character of the phase angle. Widths of the phase-angle distributions for the P wave (the time window 21 μs) are narrow for the low frequencies under 1 MHz, but they become wider in the PP P phase. This suggests that more waveform  distortions appear for longer travel distances because of reflections at the surfaces of the random medium.
Another evidence of waveform fluctuation due to random heterogeneity can be found by changes of the waveform with increasing the station distance. We refer to this type of fluctuation as a spatial fluctuation. same for the changes of station distance (station skip number), but the width of the phase distribution shows more remarkable changes with increasing the station distance.
To confirm whether or not the waveform distortion described above is actually due to random heterogeneity, we measured waveforms in steel and calculated the cross spectrum. Steel can be considered as a completely homogeneous medium. Fig. 5 shows waveforms in steel. The block size and array configuration are the same as those of Westerly granite, and the source signal is 0.5 MHz. The direct P phase appears at about 16 μs and the reflected PP P , at about 43 μs. The reflected PP P P P appears at about 71 μs. The weak waves around 30 μs correspond to the direct S phase. Waveform appears to show correlations over all observation points of the array. Fig. 6(a) shows amplitude and phase of cross spectra at the direct P wave (18 μs) and the reflected P wave (42 μs): PP P . Fig. 6(b) shows amplitude and phase of cross spectra at PP P (42 μs) for different station distances. There are strong 0.5-MHz peak and weak 2-MHz peak in amplitude spectra. The 2-MHz peak corresponds to the characteristic frequency of the PZT transducer. This indicates that we can detect waves up to 2 MHz if the medium is homogeneous. However, distribution of phase spectra becomes wider with increasing frequency, suggesting a systematic error caused by inaccuracy of beam positioning.
All the phase plots show that phase distributes in a narrow range up to 1 MHz. Figs 6(b)-1 B, (b)-2 B and (b)-3 B suggest no spatial fluctuations in waves propagating through steel up to 1 MHz. If we consider that the phase fluctuation over 1 MHz in steel is due to the errors of spectrum estimation and the errors associated with experiments (mainly beam positioning), we can conclude that reliability of the present experiments and analysis can be assured up to 1 MHz. Phase fluctuations in Westerly are remarkably lager than those in steel for the P and PP P phases. Spatial dependence of phase fluctuation is also more apparent in Westerly granite. These suggest that phase fluctuations in rock samples are caused by the random heterogeneity inherent to the rock.  between −π and π. Dotted lines indicate the upper and the lower boundaries within which 2/3 of the whole data appear in both sides of the median value. The width of the boundary roughly corresponds to the standard deviation of the Gaussian distribution covering about 68 per cent of the whole data. We can use the width of the 2/3-data boundaries to characterize the phase distribution.

Distribution of phase spectra
In Figs 7(a) and (b), histograms of the data are also shown. Symbols in the histograms indicate the 2/3-data boundaries. The phase distribution for 0.5-MHz becomes narrow when the signal waves, P, S and PP P appear. However, the phase distribution for 1 MHz becomes wide for all the lapse time except for the first P-phase. This indicates that the scattered waves are strong in 1 MHz and the signals of S and PP P are masked by the scattered waves. The distributions seem close to the Gaussian distribution when the widths of 2/3-data boundaries are small.
Due to a cyclic character of phase angle, the Gaussian probability density function (PDF) over (−∞, ∞) can be converted to the wrapped normal distribution determined in a range (−π , π): where σ 2 is the variance of the unwrapped normal distribution. Eq. (15) means that the tails of the Gaussian PDF are folded at π ± 2nπ , (n = 0, 1, 2, . . . ) and added within the range (−π, π ). Eq. (15) is shown in Fig. 8 for different σ values. When σ is small, PDF quickly decreases to zero, and the most of phase values appear  in (−π, π). However, when σ becomes large, PDF becomes flat in (−π , π), showing that PDF is close to the uniform distribution. The distribution becomes almost uniform when σ = π. Fig. 9 shows the upper and the lower boundaries of the 2/3-data region as a function of frequency at the time window 58 μs for the waves generated from the 1-MHz source. This time window corresponds to the arrival of PP P phase. The boundaries of the 2/3-data region expand as frequency and the station distance increase. The boundaries are almost the same when the station distance becomes 5 or larger. The width of the distribution is narrow below 0.75 MHz, but it suddenly becomes wide as frequency increases higher than 0.75 MHz. Then PDF of the phase distribution becomes almost uniform over 0.75 MHz when the station distance is larger than 5.
Figs 10(a), (b) and (c) show the 2/3-boundaries of phase fluctuation for the pairs with the station distance 7 as a function of time at 21, 24, 44 and 58 μs for the waveforms generated by 0.25-, 0.5-, 1-and 2-MHz source signals. The time windows at 21, 24, 44 and 58 μs correspond to the following waves, respectively: just after the direct P arrival, the tail of P wave, the scattered wave part without signal waves, and the reflected PP P . At 21 μs, phase fluctuations are small below 1 MHz, suggesting the small fluctuations in the direct P phases. At 24 μs, phase fluctuations become large at around 0.5-1 MHz depending on the source-signal frequency. This suggests that the tail of the direct P phase contains fairly large amount of scattered waves. At 44 μs, the widths of the 2/3-data range are large, showing weak correlations between waves. At 58 μs, phase fluctuations are almost similar to that at 24 μs, suggesting the signals of PP P phase.

Effect of source radiation pattern for different signal frequencies
Phase fluctuation differs for each source-signal frequency. The difference may be caused by the different radiation patterns for different source-signal frequencies. Source radiation pattern depends on the Poisson's ratio of the medium and the source size to wavelength ratio (Tang et al. 1994). Fig. 11 shows P-and S-wave radiation patterns for the frequencies 0.25, 0.5, 1 and 2 MHz. In the figures, radiation patterns are plotted on a plane which includes the rotational symmetric axis since the patters are axisymmetric. The radiation pattern was calculated for a disc-shaped compression transducer by using equations given by Tang et al. (1994), where a harmonic vibration of the vertical stress is uniformly distributed inside a circular area (a piston force). The figures show the intensity of the P and S waves from a disc-type source. The radiation pattern has maximum P-wave energy along the axial direction of the disc. The P-wave patterns show lobes having more round shape for low frequencies; 0.25 Figure 11. Energy radiation pattern of the present PZT source calculated from Tang et al. (1994). Strong S waves appear when signal frequencies are 0.25 and 0.5 MHz. This may cause strong scattered waves in the waveform around PP P for the source signal 0.5-MHz. and 0.5 MHz. The S-wave radiation pattern has a node in the axial direction, and its intensity is stronger in low frequencies; 0.25 and 0.5 MHz. Strength of scattered waves depends on scattering ability of the medium and the energy radiation pattern from the source. Since scattered wave energy is generally larger in the lateral and the backward directions (Sato 1984;Sato & Fehler 1998), the stronger radiation of P and S waves in the lateral direction causes stronger scattered wave in the medium. However, strength of the scattered wave also depends on the wave frequency: high-frequency waves are more scattered than the low-frequency waves.
For a 0.25-MHz source, energy radiation in the lateral direction is large but scattering is not strong because the wavelength is larger than the heterogeneity scale a. However, for 0.5-MHz source, energy in lateral directions is still large and scattering is stronger than the 0.25-MHz case. For 1-and 2-MHz sources, scattering becomes even stronger, but laterally radiated energy becomes smaller, and most of the scattered waves are generated from forward scattering rather than from backward scattering. To interpret Figs 10(a)-(c) for different source frequencies, we have to take those conditions into account. The different patterns of phase fluctuations for different source frequencies may be caused by the different conditions of source radiation pattern and scattering.
Although the conditions of wave fluctuation are different between the four different-frequency experiments, the onsets of the extension of 2/3-data boundaries are clear between 0.5 MHz and 1.5 MHz. This suggests that the phase angle in the seismic signal is distorted depending on the wave frequency, and the distortion suddenly becomes strong between 0.5 and 1.5 MHz.

Phase angle fluctuation against distance
Waves received at closer stations are more correlate each other, even if they have fluctuations due to scattering. Waves become less correlated as station distances increase. Correlation of waves depends on the wavelength of the seismic wave and the scale length of random heterogeneity. Since the width of the 2/3-data region is associated with the shape of the phase distribution function, the width indicates incoherency of the seismic wave. We can then describe the relation between seismic wavelength and the scale length of random heterogeneity by plotting the width as functions of frequency and the station distance of the circular array. A real distance d between stations is obtained by the station-pair number N between stations: where r is the radius of the circular array (10 mm). When plotting the phase fluctuation against distance and frequency, it is better to use the parameter normalized by the wavenumber k. The wavenumber is calculated from the seismic velocity of Westerly granite v(4.78 km s −1 ) as 2π f /v = 2π/λ, where f and λ are the frequency and the wavelength of P wave, respectively. We use the parameters ka and kd, which are the products of the wavenumber and the correlation distance a of the random medium and the space distance d, respectively. ka values are shown in Table 1. Once we describe the relationship between wave fluctuation and the random heterogeneity by using ka and kd, the relationship holds for any scales; from laboratory models to the real earth. Figs 12(a) and (b) indicates the normalized relationship for the 0.5-and 1-MHz source signals. They show a sudden increase of the width of the 2/3-data range when ka becomes 0.2-0.3, which suggests a critical frequency. The change of the 2/3-data width with respect to kd indicates the wave correlation against station distances. In Figs 12(a) and (b), low ka values show small values of the 2/3-data width, indicating that waves well correlates in lower frequencies. For the ka values larger than 0.2-0.25, the width increases with increasing kd values when kd < 5, suggesting that waves become uncorrelated with increasing station distance as shown in Fig. 9. Then the width becomes independent of kd at kd > 5, indicating the critical distance over which waves become fully incoherent. We can see the area of correlated waves in random media, in terms of space distance and seismic wave frequency. Figure 12. Width of 2/3-data range plotted on the ka−kd plane at the lapse time of PP P arrival. The width rapidly increases when ka becomes 0.2-0.3. The width depends on kd for kd < 5 but it becomes independent for kd > 5. The area of random phase fluctuation can be assigned as high-valued area of the contour plot.

Wave fluctuation in random media
In the frequency-domain reconstructions of subsurface velocity structure, spectrum estimation is important, but it becomes erroneous when the signal waves are strongly fluctuated by scattered waves. It is important to note that the fluctuated waves are not random noise in the time-series but signal-induced waves. As shown in the present analysis, amplitude fluctuation is almost independent of frequency but phase fluctuation strongly depends on frequency. Thus the relationship between phase fluctuation and frequency is more important than the amplitude fluctuation.
The strong fluctuation of waveform at the ka value about 0.2-0.3 was also pointed out by Sivaji et al. (2001b). They analysed wave data similar to that presented here; for the wave data obtained for rocks having different space correlations, a, but almost the same intensity of fluctuation, . They found that the following parameters of waveform are suddenly changed around ka = 0.2-0.3: the waveform correlation, standard deviation of traveltime fluctuation, standard deviation of the total waveform correlation, and the standard deviation of wave-energy fluctuations. Their results can be interpreted by the fact that the phase fluctuation is realized by the uniform distribution when ka value becomes close to 0.2-0.3.
To interpret experimental results by using a general relationship between seismic waves and random media, a figure describing the relationship between ka and kL is very useful, where L is the travel distance of the seismic waves. Fig. 13 illustrates a classification of scattering problems conjectured by Aki & Richards (1980) (Vol. II, p. 749), showing applicable methods in the ka−kL diagram. For drawing the figure, they tentatively assume the intensity of random heterogeneity, as √ 10 per cent. The rock used in the present experiment shows a slightly larger value for , and the area of random P Figure 13. Aki and Richards' conjecture for classifying seismic wave problems in random media drawn on the ka−kL plane (Aki & Richards 1980). The area of laboratory experiments is shown as a shaded area. ka−kL relationships of the present experiments are shown as the dashed lines corresponding to the direct P and PP P phases. media (central part) will extend downwards more. Laboratory experiments roughly cover the shaded area (Sivaji et al. 2001b). ka and kL relationships of the present study are shown by thick dashed lines and corresponding frequencies are indicated by horizontal thin dashed lines. The present study covers the border of two regions: equivalent homogeneous media and scattering media. The abrupt change around ka ∼ 0.2 − 0.3 corresponds to the boundary between the equivalent homogeneous medium and the random medium, which is characterized by a critical frequency in Fig. 13.

Effects of phase fluctuation on seismic inversion problems
When studying fault motions or underground velocity structures on the basis of seismic observations, frequency-domain inversion techniques are adopted to reconstruct fault motions or velocity structures. Langston (1979) and Ammon et al. (1990) proposed the receiver function to reconstruct underground velocity structures. Their method calculates spectra for the components of seismic particle motions and obtains response functions in the frequency domain. Pratt (1990) and Pratt et al. (1996) have applied frequency-domain waveform inversion technique for reconstructing velocity structures between seismic sources and receivers of two boreholes. Cotton & Campillo (1995a) and Cotton & Campillo (1995b) used a frequencydomain inversion algorithm to reconstruct fault motions from observed seismic waves. In all the above techniques, accuracy of spectra is crucial for obtaining reliable solutions. However, if we consider that the subsurface structures contain small-scale random heterogeneity, phase fluctuations of seismic waves are affected by scattering from the random heterogeneity. The present results suggest that there is a limit of resolution in phase spectra controlled by the relationship of the seismic wavelength and the scale length of random heterogeneity; ka ∼ 0.2 − 0.3. It suggests that there is a resolution limit when applying frequency-domain waveform inversion because incoherency of phase spectra may lead to incorrect results.
On the other hand, there is another reconstruction technique that uses only amplitudes information of seismic waves. Nakahara et al. (1999) and Nakahara et al. (2002) used waveform envelope inversion to reconstruct earthquake fault motions and to explain the generation of high-frequency seismic waves. They avoided difficulties in treating the phase fluctuations at high frequency and used only amplitude information that is independent of frequency. However, if we carefully notice the incoherency of the phase spectra in ka > 0.2 − 0.3, we would be able to use phase information in seismic waves up to the frequency corresponding to that ka value.
When comparing waveforms of observed waves and the waves calculated from an estimated model, it is natural to consider that observed waves are distorted due to scattering from the small-scale unresolved structures. The scattered waves are the signal-induced waves. Langston & Hammer (2001) pointed out that the signal induced waves are decoupled from components of seismic waves and cause more serious problems in receiver function analysis than spectral holes in the frequency-domain analysis. In subsurface image reconstructions, serious problems associated with waveform fluctuation also exist. The accuracy of spectral estimation in deconvolution techniques is similarity affected by waveform fluctuation. These problem can be addressed by noticing the distribution of phase fluctuation as a function of ka value.
Subsurface reconstruction images are truncated to a certain wavelength in space (Snieder & Trampert 1999;Spetzler & Trampert 2003), and only the large-scale heterogeneities are shown. The large-scale heterogeneities are expressed by blocks, by truncated functional expansions, or by spline-smoothed curves. Dividing into small blocks merely enhances high frequency artefacts in the resulting image corresponding to the smaller split blocks, and does not approach reality. It is thus misleading to expect that dividing into smaller cells (or merely expanding structures to higher order functions) provides a more realistic image. We must consider that seismic waveform fluctuation depends on ka value, and the misfits between calculated waves and observed waves also depend on ka. The misfits are thus related to the small-scale heterogeneity due to the scattering from the small-scale randomness of subsurface structures (Snieder & Trampert 1999;Spetzler & Trampert 2003).

C O N C L U S I O N S
We have presented experimental results about waveform fluctuation in a small-scale random medium. Since seismic waves are strongly affected by the small-scale heterogeneity (short-wavelength structures), we must be careful about selecting wave frequencies for underground imaging or for reconstructing fault motions by seismic waves. Amplitude fluctuation is almost independent of frequency but phase fluctuation is strongly dependent on frequency. At low frequencies, phase fluctuations are close to Gaussian distributions with narrow distribution width. At high frequencies, phase fluctuation is random, with an apparently uniform distribution between −π and π, causing incoherent waves in a seismogram. This is very important when we use frequency-domain reconstruction techniques of geophysical inversion problems. In the frequency domain analysis, spectra are the complex number with amplitude and phase information. The uncertainty in phase spectra is not constant over all the frequencies, but the uncertainty changes with frequency and misfits in spectra should be treated on the basis of this kind of frequencydependent error. To realize the concept of the frequency-dependent error, knowledge about the distribution of the misfits is of primary importance. The present experimental results about phase fluctuation provide a plausible distribution of phase fluctuation in random media, and help to improve frequency-domain reconstruction techniques of geophysical inversion problems.

A C K N O W L E D G M E N T S
The authors thank Haruo Sato and Tatsuhiko Saito for helpful discussions and comments. Robert Kranz also gave us valuable comments for improving the manuscript. Comments from two anonymous reviewers were also very important and helpful.