Docstoc

C seismology

Document Sample
C  seismology Powered By Docstoc
					                                                         1




Z-99 A NEW S WAVE IMAGING METHOD: A
          RECEIVER FUNCTION APPROACH IN THE
          TAU-P DOMAIN
                                                                                                P. EDME, S.C. SINGH
                Laboratoire de Geosciences Marines, Institut de Physique du Globe, 4 place Jussieu 75005 Paris, FRANCE




                                               Abstract
 Recovery of both P and S wave information is important to accurately determine sub-surface elastic
 and structural properties. We propose to use the receiver function (RF) technique to analyse multi-
 component data. The RF method is traditionally employed in global seismology to enhance converted
 waves produced by earthquakes when propagating within the upper mantle and the crust. In this study
 we present a scheme that adapts this technique to the seismic exploration domain. The calculation of
 RF is based on an iterative deconvolution of seismograms in the τ-p domain. This process allows the
 separation of upgoing the P and S wavefield, and clearly highlights the main shear wave arrivals by
 computing the impulse response of the prospected section. In addition our method does not require a
 priori information and can be used for S wave AVO study as well as S wave imaging. Numerical tests
 on synthetic wide-angle multi-component OBC data show encouraging results and demonstrate the
 potential of this new method in producing accurate images of the S wave structure with higher signal
 to noise ratio.

                                              Introduction
 In recent years, there has been an increasing interest in the acquisition and analysis of multi-
 component seismic data. These data provide informations on both the P and S wave structure of the
 subsurface. However, the conventional methods developed for P waves analysis are often used to
 interpret multi-component recordings. These approaches yielded, and will no doubt continue to yield,
 important scientific results but developing new methods is necessary. Here we propose to use receiver
 functions to analyse OBC data.

                                      Receiver function generality
 Currently receiver functions are used exclusively in global seismology and consist of analysing
 earthquake signals. Main characteristics of the technique have been explained by many authors
 (Langston, 1979; Ammon, 1991) enabling the receiver function method to be routinely applied. This
 technique has now become a standard method for mapping structure and analysing shear wave
 velocities from both the vertical and horizontal component. A classical receiver function is a time
 series and is defined as the deconvolution of the horizontal component by the vertical component of
 the data, i.e
                      Ux(t) = RF(t) * Uz(t)    or       Ux(ω) = RF(ω) . Uz(ω) ,                    (1)
  where “*” is the convolution operator and “.” is the multiplication.
  A receiver function can either be thought of as the impulse response in the temporal domain or the
 transfer function in the frequency domain (equations 1), linking the S and the P wavefields recorded
 on the Ux and Uz components in the case where the teleseismic event propagates nearly vertically
 below the receiver (figure 2a). The deconvolution removes the effect of the seismic source and the
 travel path in the mantle and emphasises the P to S converted waves.
 The RF computation can be done either in the frequency domain or in the temporal domain, and
 should in theory lead to similar results. The deconvolution in the frequency domain is the most
 commonly employed and consists of dividing the Fourier transforms of Ux(ω) by Uz(ω). To avoid
 numerical problems the near zero values of Uz(ω) are replaced with a fraction of the maximum value,
 which enables the stabilization of the division. This fraction is known as the water-level. Figure 1b
 shows the result of the deconvolution in the frequency domain using a water-level equal to 1%. More


                EAGE 66th Conference & Exhibition — Paris, France, 7 - 10 June 2003
                                                    2


recently, a new deconvolution method in the time domain has been employed for the impulse response
calculation (Ligorria and Ammon 1999).
                             RF(t)=Σi ai . δ(t-ti) .                                               (2)
 This method uses the mathematical development of Kikuchi and Kanamori (Kikuchi and Kanamori
1982) and is based on the least squares minimization of the difference between the observed horizontal
component and the predicted signal generated by the convolution of an iteratively updated signal spike
train (equation 2) and the vertical component seismogram. Results of this deconvolution is shown in
figure 1c.




Figure 1: Examples of deconvolution of the seismogram in (a) by the first second of the signal. (b):
Signal obtained in the frequency domain. (c): Signal obtained by the iterative deconvolution method.
This signal spike train explains 96% of the original signal energy in (a).

In comparison with the RF in the frequency domain, the iterative deconvolution shows several
desirable qualities. The calculated RF is a true impulse response composed of a sum of Dirac delta
functions. The process uses crosscorrelation between the signal to deconvolve and a selected source to
firstly detect the biggest arrival in the seismogram and strip it before extracting smaller phases. The
iterations are repeated until a chosen percentage of the total energy has been removed. Like this, even
after computation, the interpreter has the possibility to choose the RF as a function of the degree of
accuracy and detail required. In addition, in each iteration, the quality of the deconvolution can be
evaluated by comparing the detected phase energy before and after deconvolution (this parameter is
called the goodness of fit). For these reasons, we use iterative deconvolution in the time domain to
compute our receiver functions.

                           Application to the seismic exploration problems
To adapt the receiver function technique to the seismic reflection case, we first have to separate the
upgoing P and S wavefield. This kind of operation has been studied by Wang et al. (2002) and can be
done by a simple rotation of components in the τ-p domain (where τ is the intercept time and p the
horizontal slowness). The transformation in the τ-p domain rearranges seismic records into a ideal
structure for processing. In the t-x domain, all waves are recorded at the station with large variation of
incidence angles (figure 2b), while in the τ-p domain, incidence angles of upgoing waves depend only
on the wave type (i.e P or S arrival) and do not vary over time for a given constant horizontal slowness
(figure 2c).




Figure 2: (a): Ray geometry in global seismology (seismologists select only teleseismic events
coming from 30 to 90° angular distance to be sure of near-vertical wave propagation below the
receiver). (b): Ray geometry in seismic exploration. Note the increase in variation of incidence angles
of waves in the reflection case. In contrast, in the τ-p domain (c), incidence angles are constant for
each wave type. For a given p, the P wave frame is defined by the Urp and Utp component (radial and
tangential of the direction of propagation of P waves in the shallower layer, respectively).
                                                   3


The separation scheme is based on the analysis of the first reflected phase P1P in the τ-p domain, and
uses the deconvolution technique described above to evaluate the incidence angle φp of P waves for a
given p. Considering a1(p) as the amplitudes of the Dirac delta functions calculated by deconvolving
Ux]p1p by Uz]p1p of the selected first arrival, angle φp satisfies the following equation:
                              φp(p) = - atan(a1(p))                                                (3)
 Comparison between theoretical and estimated angles as a function of the slowness reveals a good
precision except at small p, where not enough energy is present on the Ux component, but inaccuracies
can be detected by observing the goodness of fit.
The calculated incidence angles φp are then used to rotate the components into the P wave frame
composed of the orthogonal Urp and Utp components (figure 2c).
             Urp = Ux.sinφp – Uz.cosφp     and Utp = Ux.cosφp + Uz.sinφp                           (4)
 On the Urp component all P wave energy is recorded with additionnal projection of PS wave energy.
On the Utp component, which is normal to the direction of propagation of P waves, the wavefield is
free of any compressionnal energy and consequently only the shear wavefield is present. This
separated S wavefield is then iteratively deconvolved by the zero phase wavelet determined from the
first arrival P1P on the Urp component. This change in the phase spectrum of the source enables a
better quality of the deconvolution, providing a better fit with both maximum and minimum phase
wavelet.

                                           Synthetic example
We begin with synthetic wide-angle OBC shot-gathers in the time-offset domain calculated with a
finite-difference code for a simple plane layered model, consisting of two layers over a half space and
below a water layer. Free-surface multiples are not included and the direct P wave arrival is removed.
The source signal is a Ricker wavelet with peak frequency of 10 Hz. Offsets range from 20 m to 10
km with equal interval of 20 m. Both vertical and horizontal components are transformed in the τ-p
domain and a static correction is applied. Figure 3a,b show these original seismograms. Both
components contain P and S waves energy.




 Figure 3: (a) and (b) are the original Uz and Ux components in the τ-p domain. The direct arrival is
removed except on the Uz component for precision. (c) shows theoretical incidence angles of P waves
(solid line) and calculated angles (crossed curve) . Inaccuracy occurs for small values of slowness
where the goodness of fit is less than 80%. (d) shows the isolated S wavefield on the Utp component.
(e) and (f) are the RF representing 60% and 90% of the energy of the Utp component, respectively.

Analysis of the first reflected wave allows the determination of the incidence angles of the P waves
(figure 3c). Inaccuracy occurs when incidence angles are less than 10°. However, for higher values the
               EAGE 66th Conference & Exhibition — Paris, France, 7 - 10 June 2003
                                                   4


goodness of fit can be more than 99%, which shows the accuracy of this deconvolution method and
demonstrates its applicability to analyse wide-angle data. The S wavefield obtained by the rotation is
isolated on the Utp component and is plotted in figure 3d. Deconvolution is applied next and the
calculated impulse responses (the RF) are shown in figure 3e,f, which represent respectively 60% and
90% of the original Utp signal energy.

To demonstrate the potential of our method, we apply it to a 2D velocity model with the second
interface a 18° dipping. The distance between two shots is 100 m. The same process (as described in
figure 3) is performed on 47 shots. The RF is then convolved with a Gaussian wavelet to be consistent
with data and to enable constructive summation. An elliptical moveout correction is performed and τ-p
traces are stacked for p interval from 0.003 km/s to 0.3 km/s. Figure 4 shows the results of the stack
applied on the Ox, Utp, and RF signals using the same PS velocity model.




Figure 4 : Stacking sections at different steps of the process. (a): Input velocity model. (b): Results
from the stack of the horizontal component Ux which contain both the S wavefield plus P waves. (c):
The Utp component contains only S wave arrivals. (d): The RF stacked section where the chosen RF
represent 80% of the energy of Utp signal. The RF have previously been convolved with a Gaussian
function whose width is adjusted in order to respect the vertical resolution.

                                            Conclusion
Comparing with classic stacks of the Ux components (figure 4b), which are often used to map the shear
wavefield, the receiver function treatment has clearly improved the quality of mapping S wave
structure. In a first step, remaining P waves recorded on the horizontal components are removed,
leaving isolated only S wave arrivals (figure 4c). This separation is performed by deconvolution
analysis and rotation of recorded seismograms in the τ-p domain. It shows a good performance
especially when the incidence angles increase and also when the unwanted P wavefield projected on
the Ux component becomes important. The rest of the iterative deconvolution allows us to highlight
the S arrivals by both increasing the signal to noise ratio and producing a sharp resolution of the
impedance contrasts.

                                        Acknowledgements
This work is being founded by the LITHOS consortium located in Cambridge and Paris. I thank all
members of the group for their kindness and their constructive remarks.

                                              References
Ammon, C.J., 1991, The isolation of receiver effects from teleseismic P waveforms: Bull. Seism. Soc.
Am., 81, 2501-2510.
Kikuchi, M., Kanamori, H., 1982, Inversion of complex body waves: Bull. Seism. Soc. Am., 72, 491-
506.
Langston, C.A., 1979, Structure under Mount Rainier, Washington, inferred from teleseismic body
waves: J. Geophys. Res., 84, 4749-4762.
Ligorria, J.P., Ammon C.J., 1999, Iterative deconvolution and receiver function estimation: Bull.
Seism. Soc. Am., 89, 19-36.
Wang, Y., Singh, S.C., Barton, P.J., 2002, Separation of P and SV wavefields from multi-component
seismic data in the τ-p domain: Geophys. J. Int., 151, 663-672.
                                 5




EAGE 66th Conference & Exhibition — Paris, France, 7 - 10 June 2003