Bulletin of the Seismological Society of America
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Bulletin of the Seismological Society of America; February 2006; v. 96; no. 1; p. 1-10; DOI: 10.1785/0120050036
© 2006 Seismological Society of America
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF)
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Nowack, R. L.
Right arrow Articles by Sheng, J.-M.
Right arrow Search for Related Content
GeoRef
Right arrow GeoRef Citation

Article

Correlation Migration Using Gaussian Beams of Scattered Teleseismic Body Waves

Robert L. Nowack1, Saptarshi Dasgupta1, Gerard T. Schuster2 and Jian-Ming Sheng2

1 Department of Earth and Atmospheric Sciences
550 Stadium Mall Drive
Purdue University
West Lafayette, Indiana 47906
 (R.L.N., S.D.)

2 Geology and Geophysics Department
University of Utah
Salt Lake City, Utah 84112
 (G.T.S., J.M.S.)


    Abstract
 Top
 Abstract
 Introduction
 Description of the Method
 Applications of Gaussian Beam...
 Conclusions
 
Correlation migration for structural imaging using Gaussian beams is described for the inversion of passively recorded teleseismic waves. Gaussian beam migration is based on an overcomplete frame-based representation of the seismic wave field and uses localized slant-stack windows of the data. Paraxial Gaussian beams are then utilized for the backpropagation of the seismic waves into the medium. The method can provide stable imaging of seismic data in smoothly varying background media where caustics and triplicated arrivals exist. We develop Gaussian beam migration for structural imaging, which allows for incident teleseismic waves from beneath the structure, using directly scattered or surface-reflected phases. The method is applied to synthetic data computed for a collisional-zone model, and the inversions show that the Gaussian beam migration can image structure using passive teleseismic data. The method is next applied to synthetic data that have been convolved with an observed pulse from the 1993 Cascadia experiment. The autocorrelation is computed and the autocorrelated data are migrated by using Gausian beams for direct P-to-S scattered waves and surface-reflected waves. Although the migrations of the autocorrelation data with the source pulse included have more noise than the migrations of the impulsive source data, the subduction zone structure can still be clearly identified in the migrated results without removal of the source pulses.


    Introduction
 Top
 Abstract
 Introduction
 Description of the Method
 Applications of Gaussian Beam...
 Conclusions
 
Teleseismic data are now being used to image crustal and deeper upper-mantle structure using broadband seismic data from distant earthquakes. One approach is to use scattered S waves from incident teleseismic P waves to image the velocity structure beneath a seismic array. Receiver functions are constructed by using the vertical P-wave component to deconvolve the horizontal components (Vinnik, 1977; Langston, 1977, 1979; Owens et al., 1984). More recently arrays of broadband stations have been deployed using teleseismic sources to image crustal and upper-mantle structure. To image all the components of the seismic data more general processing steps have been applied, such as array stacking for the source time function (Li and Nabelek, 1999; Langston and Hammer, 2001), eigenimage analysis (Rondenay et al., 2001; Ulrych et al., 1999), auto- and cross- correlation techniques (Sheng et al., 2001, 2003; Schuster et al., 2003, 2004; Yu et al., 2003) and multichannel deconvolution (Bostock, 2004).

Initially, basic stacking techniques were used for the imaging of receiver function array data (Dueker and Sheehan, 1997; Zhu, 2000). More recently migration techniques have been applied to P-to-S conversions from teleseismic waves recorded by passive arrays (Bostock and Rondenay, 1999; Ryberg and Weber, 2000; Sheehan et al., 2000; Poppeliers and Pavlis, 2002; Pavlis, 2003). This has been extended to the imaging of surface-reflected phases by Bostock et al. (2001) in a formal ray/Born migration approach and applied to data from the Casc93 experiment by Shragge et al. (2001) and Rondenay et al. (2001). Ghost reflections from the free surface have been described by Sheng et al. (2001, 2003), Yu et al. (2003), and Schuster et al. (2003, 2004) for both exploration and teleseismic imaging and have also been used for daylight interferometric imaging with unknown or random source signals (Rickett and Claerbout, 1999). In this article, the Gaussian beam migration method is investigated for the inversion of teleseismic data with application to passive imaging experiments.

Gaussian beam migration uses an overcomplete frame of smoothly localized Gasussian windows to represent the seismic data and paraxial Gaussian beams to propagate the data back into the medium. Because, for Gaussian beam migration, the wave fields are decomposed into Gaussian beams, the imaging condition then uses individual Gaussian beams and allows for caustics as well as triplicated seismic wave fields in the background medium. In contrast, Kirchhoff migration or ray/Born inversions require the first arrivals or most energetic arrivals at the scatterer and this results in an incomplete imaging condition unless further analysis is performed. This is an advantage of Gaussian beam migration over other high-frequency migration approaches and is one of the reasons that Gaussian beam migration has become one of the principal migration tools in the petroleum industry. Nonetheless, even in the oil industry, Gaussian beam migration is primarily used for structural imaging and true amplitude formulas are still being developed (N. R. Hill, personal comm., 2003; Albertin et al., 2004).

We test the Gaussian beam approach for structural imaging using synthetic teleseismic data for a collisional zone model similar to that used by Schragge et al. (2001). In addition to testing Gaussian beam migration algorithm for passive imaging, we apply Gaussian beam migration to autocorrelation data. The synthetic data with an impulsive source pulse are convolved with an observed pulse from the 1993 Cascadia experiment data, and then migration is applied to the filtered autocorrelation data including the source pulse. Although the imaging results have more noise than for the impulsive data, clear images of the structure are still obtained. The correlation approach using Gaussian beams has the potential of using all components of the seismic data in addition to just the horizontal components for the imaging. Although true amplitude imaging has not been used for this study, based on synthetic tests in complicated background media using seismic reflection data, the Gaussian beam method can already provide improved structure images compared with existing ray/Born and Kirchhoff imaging techniques (Nowack et al., 2003).


    Description of the Method
 Top
 Abstract
 Introduction
 Description of the Method
 Applications of Gaussian Beam...
 Conclusions
 
For an incident, teleseismic plane wave from beneath the structure, the isotropic, elastic equation can be written


Formula 001

(1)
where ui is the particle displacement, {lambda} and µ are the elastic constants, and {rho} is the density (Aki and Richards, 1980). The commas signify differentiation. Assuming a constant density, and perturbing the elastic parameters as {lambda}(x) = {lambda}0(x) + {delta}{lambda}(x) and µ(x) = µ0(x) + {delta}µ(x), this results in a perturbed, particle displacement Formula . To first order, this can be written as L0[{delta}u] = –{delta}L[u0], where {delta}L is the first-order perturbation from the initial differential operator L0 with


Formula 002

(2)
where


Formula 003

Assuming further that the medium is a Poisson solid, then {delta}{lambda} = {delta}µ = {rho}{delta}β2, with {alpha} and β being the P- and S-wave speeds. The Qi term can then be expanded as


Formula 004

(3)
This is similar to equation 13.22 of Aki and Richards (1980) for a Poisson solid with constant density. For an incident planar P wave, Formula can be written at high frequencies as Formula , where T(x) is the travel time, A(x) is the amplitude, pi = {partial}T/{partial}xi is the ray parameter vector perpendicular to the incident teleseismic wavefront and Formula . The Qi term resulting from scattering is then approximately


Formula 005

(4)

In the frequency, domain equation (2) with the body force term Qi from equation (4) can be rewritten as


Formula 006

(5)
where the integration is over the subsurface scattering locations. The Green's function is gin (x', xg, {omega}) with the locations of the sources placed at the geophone locations xg and with source components n in the background medium. In the far field, the approximate high-frequency Green's function in the background medium can be written


Formula 007

(6)
where this includes the individual P and S phases propagating downward from the array to the scatterers. The incident wave Formula used for imaging here will be either an upward- propagating direct P wave or a surface-reflected P and then scattered as a P or S wave.

Approximating equation (5) with the highest-order frequency term for Qi in equation (4) and using equation (6) results in


Formula 008

(7)
where the dependence on the incident horizontal wavenumber components ps of the incident source wave field has been included. Also, the source time function S({omega}) has been incorporated with Formula signifying the incident wave field for an impulsive source. The terms with the derivatives of {delta} µ can be written in a similar fashion, but for the purposes of structural imaging, I will consider only the first highest frequency term. Alternatively, Cerveny (2000) incorporated the derivative terms by using integration by parts applied to equation (5). Using equation (7) to model the P waves on the x3 component and the S waves on the x1 component, then


Formula 009

(8a)


Formula 010

(8b)

Writing equations (8a) and (8b) as {delta}u = A{delta}µ, then adjoint operators to these can be obtained from the definition of the adjoint


Formula 011

where the parentheses indicate the inner products with respect to the variables indicated by the subscripts, the overbar indicates a scalar complex conjugate, and A* is the adjoint to A. In addition, this relation holds for homogeneous boundary conditions. The adjoint operator can be obtained from this relation in a straightforward manner as {delta}µ = A*{delta}u. For equations (8a) and (8b) the adjoints can be written


Formula 012

(9a)


Formula 013

(9b)
where {delta}ui(xg, ps, {omega}) are the observed seismic data, gi1 and gi3 are the Green's function components from the geophone positions xg to the subsurface scatterer positions x', the source time function is S({omega}), and {omega} is frequency. The overbars again indicate scalar complex conjugates. The impulsive source wave field is represented by Formula (x', ps, {omega}) for an incident plane wave from below specified by the horizontal wavenumber components ps. The complex conjugated terms can be understood as the cross-correlation of the backpropagated observed seismic data with the incident source field and is equivalent to the imaging condition given by Claerbout (1976). Equations (9a) and (9b) can then be used for migration of teleseismic data for structural imaging in either 2D or 3D. Correlation migration is a kinematic formulation of the imaging problem and can be related to the adjoint operator of equation (8). Even in the petroleum industry with large amounts of data this provides excellent structural images (Hill, 2001). However, for the imaging of passively recorded teleseismic data, the data quality may not warrant the use of more complete imaging formulas based on the data coverage currently available.

In the Gaussian beam migration approach, the Green's functions are decomposed into Gaussian beams. The summation of Gaussian beams for forward modeling is an asymptotic method for the computation of wave fields in smoothly varying inhomogeneous media, which describes high-frequency seismic wave fields by the summation of paraxial Gaussian beams (Popov, 1982; Cerveny et al., 1982; Nowack and Aki, 1984). Reviews have been given by Cerveny (1985a, 1985b), Babich and Popov (1989), and more recently by Nowack (2003). An advantage of the summation of Gaussian beams for constructing more general wave fields is that the individual Gaussian beams have no singularities along their paths. Also, no two-point ray tracing is required and triplicated arrivals are naturally incorporated into either the forward or inverse modeling.

A criticism of early implementations of Gaussian beam summation was given by White et al. (1987) in terms of completeness and accuracy of the summations. More recently, however, overcomplete frame-based summations of Gaussian beams have been developed based on windowed and wavelet transforms to address issues of completeness. With Gaussian windows, this was referred to as frame-based Gaussian beam summation by Lugara et al. (2003). In an overcomplete frame-based approach, the wave field is decomposed into beam fields that are localized both in position and direction (Fig. 1a) and then propagated into the medium. In this type of decomposition, position plays the role of time and wavenumber plays the role of frequency in a time- frequency style decomposition. Gabor originally suggested using modulated and translated Gaussian windows for windowed Fourier analysis. Although a basis cannot be formed by using a windowed Gabor frame, an overcomplete frame expansion can be constructed that has the additional benefit of providing redundancy in the expansion (Feichtinger and Strohmer, 1998). Overcomplete representations of seismic wave fields in terms of coherent states and related Gabor windowed transforms were also described by Thomson (2001), and Thomson (2004) investigated the seismic head wave problem by an overcomplete Gabor representation of beams.


Figure 001
View larger version (23K):
[in this window]
[in a new window]

 
Figure 1. (A) Frame-based Gaussian beam decomposition. Initial wave fields are decomposed into Gaussian windowed beam components that are then launched at different angles into the medium. (B) Example of the backpropagation of a single Gaussian beam with a low-frequency initial beam width of 30 km for a pPp phase for the model shown in Figure 2.

 


Figure 002
View larger version (15K):
[in this window]
[in a new window]

 
Figure 2. A collisional zone model.

 
An early framelike decomposition of seismic wave fields without using the language of frame theory was given by Hill (1990, 2001). By physical reasoning Hill derived an overcomplete discrete sampling of the seismic wave field and the component beams were then backpropagated into the subsurface using paraxial Gaussian beams. The backpropagation of surface seismic data from a single Gaussian beam for P scattering of an incident teleseismic P wave is shown in Figure 1b for the model shown in Figure 2. The initial beam width is approximately 30 km. In a related effort, Wu developed beamlet migration using Gabor- Daubechies frame operators to decompose observed seismic data and propagation using complex phase screens (see, for example, Wu and Chen, 2001, 2002). Although the relation between beamlets and paraxial Gaussian beams is still being investigated, the use of paraxial beams may be a faster high- frequency alternative. Recent examples of Gaussian beam migration for reflection data include Nowack et al. (2003) and Gray (2005).

The Gaussian beam expansion for the elastic Green's function in a heterogeneous medium with initial position xg and final position x' can be written (Cerveny et al., 1982; Cerveny, 1985a, 1985b; Cerveny, 2000)


Formula 014

(10)
where Formula is the individual Gaussian beam, {Psi}({gamma}, {omega}) is the weight function for the summation of the beams, and {gamma} specifies the ray coordinates for the beams. The form is similar in 2D and 3D with the differences in the weights for the individual Gaussian beams (Cerveny, 2000). This can be used to describe either the far field P- or S-wave component of the Green's function. Each paraxial Gaussian beam can be written as


Formula 015

(11)
where Formula (x', xg, {gamma}) are the Gaussian beam amplitudes including the geometric spreading, any R/T coefficients, the ray-dependent source strength and the polarization vector at the positions along the beam's central ray. The phase function T(x', xg, {gamma}) is paraxially approximated for positions off the central ray; therefore, two-point ray tracing is not required. The phase along the central ray is real, and the curvature matrix of the wavefront for a paraxially approximated Gaussian beam is complex and positive definite. The real part represents the curvature of the wavefront and the imaginary part tapers the amplitude away from the central ray forming a beam solution. The geometric spreading is also complex, as well as nonsingular and regular along the entire beam. The beam parameters are specified by the real and imaginary part of the complex curvature at some point along the beam and represent the beam width and wavefront curvature at that point. The dynamic ray equations can then be used to compute these values at other points along the ray (Cerveny, 2000). The beam parameters are commonly specified at either the initial or end point of the beam.

Following Hill (2001), I specify the initial beam parameters with planar wavefronts at array positions along the surface. If the initial point of the Green's function in equation (10) is at a slightly shifted location along the surface, a phase compensation term can be used to adjust the Green's function as:


Formula 016

(12)
where the beams are now being specified by the initial horizontal wavenumber components pg. The phase term in equation (12) performs a slight adjustment of the phase from position xL to position xg along the geophone array.

For simplicity, I will describe the 2D case for the scattered S waves on the horizontal component for an incident teleseismic plane P wave. The scattered P wave on the vertical component and the 3D extensions are analogous. The division of unity formula for Gaussian functions in the 2D case can be written as


Formula 017

for xL = L{Delta}L and {Delta}L << 2{sigma}, where {sigma} is the width of the Gaussian functions at the surface. Hill (2001) provides a similar formula for the 3D case. Inserting the division of unity formula into equation (9a), using the far-field S-wave component of the Gaussian beam expansion of the Green's function in equation (12), and assuming a single incident source wave represented by Formula results in


Formula 018

(13)
where xL is the discretely sampled beam center, and Cl({omega}) = ({omega}/{alpha}0)2 ({Delta}L)/Formula with {sigma} given below at a low- frequency reference level, and Formula are the individual Gaussian beams for the far-field scattered S-wave Green's function. The overbars again signify complex conjugates. Also,


Formula 019

(14)
are the local horizontal slant stacks of the observed data Formula in the xg coordinate that have been Gaussian windowed at the distances xL = L {Delta}L and at the dip angle related to Formula . A similar formulation is obtained for equation (9b). The terms that downward propagate the windowed slant-stack data are just the paraxial Gaussian beams Formula for the S-wave component of the Green's function initiating from the beam-center locations. The resulting formulas are similar to those given by Hill (2001) for prestack migration where here we have used an incident planar source from beneath the structure instead of a common- offset, surface source geometry.

The incident planar wave field from the source Formula could be decomposed into Gaussian beams but could also be computed by other methods. For the examples shown here we assume a teleseismic source, in which case the source can be approximated by a plane wave incident from beneath the structure with ray tracing used to obtain the approximate incident wave field at the scattering points.

An important aspect of Gaussian beam migration is that localized beam stacks at given array locations are used for the Gaussian beam migration. Thus, if the original geophone spacing is relatively dense but irregular, this could be ac counted for in the beam stack in a manner such that the beam centers xL = L{Delta}L are still regularly spaced. Issues related to trace interpolation and aliasing were investigated by Neal and Pavlis (1999, 2001) by using a stacking procedure to smooth and interpolate the wave-field data with a Gaussian smoothing window for a given slowness vector to align the traces. In the present formulation, a local beam stacking can be used to form a regular array of beam centers for the Gaus sian beam migration as well.

For a Gaussian beam specified at the initial position on the ray, the beam width at a given frequency {omega} can be written as {sigma} = W0({omega}r/{omega})1/2, where W0 is the width of the Gaussian beam at a low-reference frequency {omega}r. The spacings of the beam centers in xL and ray parameter Formula depend on the frequency of the data. In 2D, a choice given by Hill (1990) is


Formula 020

where {omega}r is the low-reference frequency for the data and {omega}H is the highest effective frequency in the data (Hill, 1990). Formula must also be adequately sampled to avoid aliasing with


Formula 021

(see Hill, 1990, equation 34b, or Hale, 1992, equation 15).

From wavelet and window analysis, for a critically sampled Gabor frame, then {Delta}L{Delta}k = 2{pi} (Feichtinger and Strohmer, 1998). However, for a stable reconstruction of the wave field recorded at the geophones, the beams must be oversampled such that {Delta}L{Delta}k < 2{pi} (Daubechies, 1992). For Formula , and Formula for a fixed {omega}, then with the previous local beam sampling in xL and Formula , this results in {Delta}L{Delta}k = ({pi}/2)({omega}/{omega}H). Thus, at the highest frequency in the data with {omega} = {omega}H, the decomposition frame used in the Gaussian beam migration algorithm is oversampled by a factor of 4 and is even more oversampled at lower frequencies. This results in a stable decomposition of the source wave field into Gaussian beams with a fixed set of paraxial beams on a location and ray parameter lattice used for all frequencies. Also, an oversampling of 4 or more results in the Gaussian functions and their dual functions being essentially the same. Even fewer beams can be used provided that the earlier inequality is satisfied, but then the dual functions must be used (Feichtinger and Strohmer, 1998). Using the values of W0, {Delta}L, and Formula , an image for each source can be numerically evaluated by the backpropagation of the beams down into the medium.

Alternative but related forms of asymptotic synthesis and inversion of seismic data include wavepath migration (Sun and Schuster, 2001), coherent states (Thomson, 2001, 2003; Albertin et al., 2001), and the Maslov method (Chapman and Drummond, 1982; Xu and Lambere, 1998). Other early applications of Gaussian beam lattice expansions for migration include Costa et al. (1989), Lazaratos and Harris (1990), and Alkhalifah (1995).

Equations (9a) and (9b) are similar to a standard diffraction stack in exploration migration. For structural imaging, Schuster et al. (2003, 2004) used stationary-phase arguments to obtain migration imaging formulas. In the following examples, we use the previous adjoint formulas incorporating Gaussian beams, which are correct in the leading-order phase terms. We term this correlation migration using Gaussian beams. For true-amplitude migration, complete weighting factors also need to be incorporated in addition to the terms in the adjoint inversion formulas. An initial effort to obtain a true-amplitude Gaussian beam migration for the reflection case was presented by Albertin et al. (2004). For the reflection case, approximate Kirchhoff weights were given by Zhang et al. (2000), and a deconvolution approach operating directly on the adjoint image was presented by Hu and Schuster (2000). Nonetheless, Hill (1990, 2001) and Nowack et al. (2003) showed that Gaussian beam migration images for structural imaging can still provide excellent images compared with Kirchhoff migration images.


    Applications of Gaussian Beam Migration
 Top
 Abstract
 Introduction
 Description of the Method
 Applications of Gaussian Beam...
 Conclusions
 
To illustrate the Gaussian beam algorithm for passively recorded teleseismic data, a synthetic model has been derived based on the model of Shragge et al. (2001) for an idealized collisional zone model shown in Figure 2. Vertical and radial component ray/Born synthetics have been generated and are shown in Figure 3, where the scattering is incorporated by using point scatterers within a background two-layer model similar to that used by Shragge et al. (2001) for their test migrations. More complete finite-difference simulations have been performed by Nowack et al. (2004), and similar inversion results were found for the Gaussian beam imaging. The source is a plane P-wave incident from beneath the structure at a 20° angle from the vertical with a Gaussian pulse width of 1 sec. The seismograms in Figure 3 have been laterally shifted by using the incident angle and are displayed in a form similar to those of Shragge et al. (2001). The modeled phases include the direct ps phase and the surface-reflected pPp, pPs, pSp, and pSs phases. Note that for simplicity the initial P has been left off of the designations for the surface-reflected phases. Thus, the pPp indicates the wave that travels to the free surface as a P wave, reflects at the free surface as a P, and gets scattered back as a P. In Figure 3, the direct P arrival will arrive at zero time but has been muted out.


Figure 003
View larger version (58K):
[in this window]
[in a new window]

 
Figure 3. Ray-Born synthetics for the collisional zone model in Figure 1. The source P wave is incident at 20° from the vertical. The scattering is incorporated by using point scatterers embedded in a two-layer velocity structure. The modeled phases include the direct ps, and the surface-reflected pPs, pSp, pSs, and pPp phases. The seismograms have been laterally shifted to correct for the incident angle. The direct P phase would arrive at time zero but has been muted out.

 

In Figure 4, Gaussian beam migration has been applied to different phase types in Figure 3. For these examples a beamwidth W0 of 30 km was used for the low-reference frequency {omega}r of 0.1 Hz, and a high-frequency {omega}H of 2 Hz was also specified. Nonetheless, a range of values of W0 could be used following the relationships given by Hill (1990) by using physical arguments and the general frame inequalities given previously. In Figure 4a, the surface- reflected pPp phase on the vertical component has been imaged and the subduction zone and Moho can be seen. The small artifacts are from other phases included in the synthetic data not imaged by this particular imaging condition. Figure 4b shows the image of the directly scattered Ps phase on the radial component and is also well imaged. The pPs surface- reflected phase is migrated in Figure 4c, and in this case the nonimaged ps phase is still visible at the shallower depths. The pPs image in Figure 4c is less stretched vertically than the image of the ps phase in Figure 4b and therefore provides better vertical resolution. To reduce the artifacts in the individual images in Figure 4, the images have been multiplicatively combined; the result is shown in Figure 4d. However, a more comprehensive procedure based on combining the images by a coherence-weighted stack has been developed by Sheng et al. (2003).


Figure 004
View larger version (68K):
[in this window]
[in a new window]

 
Figure 4. Gaussian beam migration of synthetic receiver functions. (A) Gaussian beam (GB) migration of the surface-reflected pPp phase on the Z component. (B) GB migration of the ps-direct P-to-S scattered phase on the X component. (C) This shows the GB migration of surface reflected pPs phase on the X component. Note the better depth resolution compared with the migration of the ps phase. (D) Product of the ps and pPs GB migrations.

 

One of the problems in the imaging of teleseismic wave fields is the separation of the source wave field from the scattered wave field. As an example of including a source wave field, an observed trace from the 1993 Cascadia experiment was convolved with the synthetic data shown in Figure 3. The resulting seismograms are shown in Figure 5 and for this case are not strictly minimum phase. Figure 5a displays the horizontal component and Figure 5b displays the vertical component of the synthetic data convolved with a source pulse. Figure 5c and d shows the results of the filtered autocorrelation of the traces in Figure 5a and b. In these plots only the positive lags are shown and the zero-lag signals have been muted out. As a result of incorporating the primary wave, the autocorrelation has the effect of duplicating the scattered waves while also compressing the source time function.


Figure 005
View larger version (99K):
[in this window]
[in a new window]

 
Figure 5. Synthetic seismic data convolved with an observed seismic pulse from the 1993 Cascadia experiment. In these plots, the synthetic data, including the direct wave, are convolved with an observed vertical component trace from the 1993 Cascadia data set. (A) The X component record section is shown convolved with the observed trace. (B) This shows the Z component record section convolved with the observed trace. (C) This shows the filtered autocorrelation of the convolved X component. (D) This shows the filtered autocorrelation of the Z component.

 

Figure 6 gives the results of a Gaussian beam migration of the filtered autocorrelation data in Figure 5. Figure 6a shows the migration results of imaging the pPp phase on the vertical component and Figure 6b shows the migration results of imaging the pPs phase on the horizontal component. In both cases, the images of the subduction zone structure can be easily seen but are not as clear as in the impulsive source case presented earlier. This is expected, because for this case no prior source separation was performed. Although source field separation techniques such as multichannel deconvolution can also be performed, seismic imaging and migration methods are still required on the processed seismic data to construct a subsurface image.


Figure 006
View larger version (49K):
[in this window]
[in a new window]

 
Figure 6. Gaussian beam migrations of the Z component and X component filtered autocorrelations. (A) Gaussian beam migration of the Z component autocorrelation for the pPp phase. The result is lower frequency than the equivalent receiver function migration, but the primary parts of the structure can be seen. (B) Gaussian beam migration of the X component of the autocorrelation for the pPs phase.

 


    Conclusions
 Top
 Abstract
 Introduction
 Description of the Method
 Applications of Gaussian Beam...
 Conclusions
 
In this study, correlation migration using Gaussian beams has been described for the imaging of passively recorded teleseismic waves. Gaussian beam migration is based on an overcomplete frame-based representation of the seismic wave field that uses paraxial Gaussian beams for the backpropagation of the seismic waves and utilizes localized slant-stack windows of the data. This can provide stable imaging of seismic data in smoothly varying background media where caustics and triplicated arrivals can exist. We developed Gaussian beam migration for structural imaging that allows for teleseismic waves incident from beneath the structure to image directly scattered or surface-reflected phases. We applied the method to synthetic data computed for a collisional-zone model and the results show that the Gaussian beam migration can image structure using passively recorded teleseismic waves. We then applied the Gaussian beam imaging method to synthetic data that were convolved with an observed pulse from the 1993 Cascadia experiment. The filtered autocorrelation is computed and then imaged using Gaussian beam migration for direct P-to-S scattered waves and surface-reflected waves. Although the migrations of the autocorrelation data with the source time function included have more noise than the migrations of the impulsive source data, the subduction zone structure can still be clearly identified in the migrated autocorrelation results.

Manuscript received March 4, 2005

Aki, K., and P. Richards (1980). Quantitative Seismology, W. H. Freeman, San Francisco.

Albertin, U., D. Yingst, and H. Jaramillo (2001). Comparing common- offset Maslov, Gaussian Beam, and Coherent state migration, in 71st Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 913– 916.

Albertin, U., D. Yingst, and P. Kitchenside (2004). True-amplitude beam migration, in 74th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts,949 –952.

Alkhalifah, T. (1995). Gaussian beam depth migration for anisotropic media, Geophysics60 ,1474 –1484.[CrossRef][ISI][GeoRef]

Babich, V. M., and M. M. Popov (1989). Gaussian summation method (Review), Izv. Vyssh. Uchebn. Zaved. Radiofiz. 32,1447 –1466 (translated in Radiophys. Quantum Electron. 32,1063 –1081, 1990).[CrossRef]

Bostock, M. G. (2004). Green's functions, source signatures, and the normalization of teleseismic wave fields, J. Geophys. Res., 109,B03303 .[CrossRef]

Bostock, M. G., and S. Rondenay (1999). Migration of scattered teleseismic body waves, Geophys. J. Int.137 ,732 –746.[CrossRef]

Bostock, M. G., S. Rondenay, and J. Shragge (2001). Multiparameter two- dimensional inversion of scattered teleseismic body waves. 1. Theory for oblique incidence, J. Geophys. Res.106 ,30,771 –30,782.[CrossRef]

Cerveny, V. (1985a). Gaussian beam synthetic seismograms, J. Geophys.58 ,44 –72.

Cerveny, V. (1985b). The applications of ray tracing to the numerical modeling of seismic wave fields in complex structures, in Seismic Shear Waves, G. Dohr (Editor), Geophysical Press, London, 1–124.

Cerveny, V. (2000). Seismic Ray Theory, Cambridge University Press, New York.

Cerveny, V., M. M. Popov, and I. Psencík (1982). Computation of wavefields in inhomogeneous media—Gaussian beam approach, Geophys. J. R. Astr. Soc. 70,109 –128.

Chapman, C. H., and R. Drummond (1982). Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory, Bull. Seism. Soc. Am.72 ,S277 –S317.[Abstract/Free Full Text]

Claerbout, J. F. (1976). Fundamentals of Geophysical Data Processing: With Application to Petroleum Prospecting, McGraw-Hill, New York.

Costa, C. A., S. Raz, and D. Kosloff (1989). Gaussian beam migration, in 59th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts,1169 –1171.

Daubechies, I. (1992). Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania.

Dueker, K. G., and A. F. Sheehan (1997). Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track, J. Geophys. Res.102 ,8313 –8327.[CrossRef][ISI]

FeichtingerH. G. and T. Strohmer (Editors) (1998). Gabor Analysis and Algorithms: Theory and Applications, Birkhauser, Boston.

Gray, S. H. (2005). Gaussian beam migration of common-shot records, Geophysics70 ,S71 –S77.[CrossRef][ISI][GeoRef]

Hale, D. (1992). Migration by the Kirchhoff, slant stack and Gaussian beam methods, CWD-121, Center for Wave Phenomena, Colorado School of Mines, Golden, Colorado.

Hill, N. R. (1990). Gaussian beam migration, Geophysics 55,1416 –1428.[CrossRef][ISI][GeoRef]

Hill, N. R. (2001). Prestack Gaussian beam depth migration, Geophysics66 ,1240 –1250.[CrossRef][ISI][GeoRef]

Hu, J., and G. T. Schuster (2001). Prestack migration deconvolution, in 71st Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, Tulsa, Oklahoma,984 –987.

Langston, C. A. (1977). Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves, Bull. Seism. Soc. Am. 67,713 –724.[Abstract/Free Full Text]

Langston, C. A. (1979). Structure under Mount Rainier, inferred from teleseismic body waves, J. Geophys. Res.84 ,4749 –4762.[CrossRef][ISI]

Langston, C. A., and J. K. Hammer (2001). The vertical component P-wave receiver function, Bull. Seism. Soc. Am. 91,1805 –1819.[Abstract/Free Full Text]

Lazaratos, S. K., and J. M. Harris (1990). Radon transform/Gaussian beam migration, in 60th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts,1314 –1317.

Li, X. Q., and J. L. Nabelek (1999). Deconvolution of teleseismic body waves for enhancing structure beneath a seismic array, Bull. Seism. Soc. Am.89 ,190 –201.[Abstract/Free Full Text]

Lugara, D., C. Letrou, A. Shlivinski, E. Heyman, and A. Boag (2003). Frame-based Gaussian beam summation method: theory and applications, Radio Sci.38 , no. 2,8026 .[CrossRef]

Neal, S. L., and G. L. Pavlis (1999). Imaging P to S conversions with multichannel receiver functions, Geophys. Res. Lett. 26,2581 –2584.[CrossRef][ISI][GeoRef]

Neal, S. L., and G. L. Pavlis (2001). Imaging P-to-S conversions with broadband seismic arrays using multichannel, time-domain deconvolution, Geophys. J. Int.147 ,57 –67.[CrossRef]

Nowack, R. L. (2003). Calculation of synthetic seismograms with Gaussian Beams, Pure Appl. Geophys.160 ,487 –507.[CrossRef]

Nowack, R. L., and K. Aki (1984). The two-dimensional Gaussian beam synthetic method: testing and application, J. Geophys. 89,7797 –7819.

Nowack, R. L., W. P. Chen, U. Kruse, and S. Dasgupta (2004). Imaging offsets in the Moho: synthetic tests using Gaussian beams with teleseismic waves, EOS Trans. AGU85 ,F1359 .

Nowack, R. L., M. K. Sen, and P. L. Stoffa (2003). Gaussian beam migration for sparse common-shot and common-reciever data, in SEG Int. Exposition and 73rd Annual Meeting,1114 –1117.

Owens, T. J., G. Zandt, and S. R. Taylor (1984). Seismic evidence from an ancient rift beneath the Cumberland Plateau, Tennessee: A detailed analysis of broadband teleseismic P waveforms, J. Geophys. Res. 89,7783 –7795.[CrossRef][ISI]

Pavlis, G. L. (2003). Imaging the earth with passive seismic arrays, Leading Edge22 ,224 –231.[Free Full Text]

Popov, M. M. (1982). A new method of computation of wave fields using Gaussian beams, Wave Motion4 ,85 –97.[CrossRef][ISI]

Poppeliers, C., and G. L. Pavlis (2003). Three-dimensional, prestack, plane wave migration of teleseismic P-to-S converted phases: 1. Theory, J. Geophys. Res.108 , 2112, doi10.1029/2001JB000216 .[CrossRef]

Rickett, J., and J. Claerbout (1999). Acoustic daylight imaging via spectral factorization: helioseismology and reservoir monitoring, Leading Edge,957 –960.

Rondenay, S., M. G. Bostock, and J. Shragge (2001). Multiparameter two- dimensional inversion of scattered teleseismic waves. 3. Application to the Cascadia 1993 data set, J. Geophys. Res. 106,30,795 –30,807.[CrossRef]

Ryberg, T., and M. Weber (2000). Receiver function arrays: a reflection seismic approach, Geophys. J. Int.141 ,1 –11.[CrossRef]

Schuster, G. T., F. Followill, F. Katz, J. Yu, and Z. Liu (2003) Autocorrelation migration: theory, Geophysics 68,1685 –1694.[CrossRef][ISI][GeoRef]

Schuster, G. T., J. Yu, J. Sheng, and J. Rickett (2004). Interferometric/ daylight seismic imaging, Geophys. J. Int. 157,838 –852.[CrossRef]

Sheehan, A. F., P. M. Shearer, H. J. Gilbert, and K. G. Dueker (2000). Seismic migration processing of P-SV converted phases for mantle discontinuity structure beneath the Snake River Plain, Western United States, J. Geophys. Res.105 ,19,055 –19,065.[CrossRef]

Sheng, J., G. T. Schuster, and R. L. Nowack (2001). Imaging of crustal layers by teleseismic ghosts, EOS Trans. AGU 82.

Sheng, J., G. T. Schuster, J. C. Pechman, K. L. Pankow, and R. L. Nowack (2003). Coherence-weighted wavepath migration for teleseismic waves, EOS Trans. AGU84 ,F5991 .

Shragge, J., M. G. Bostock, and S. Rondenay (2001). Multidimensional inversion of scattered teleseismic body waves. 2. Numerical examples, J. Geophys. Res.106 ,30,783 –30,793.[CrossRef]

Sun, H., and G. Schuster (2001). 2-D wavepath migration, Geophysics 66,1528 –1537.[CrossRef][ISI][GeoRef]

Thomson, C. J. (2001). Seismic coherent states and ray geometrical spreading, Geophys. J. Int.144 ,320 –342.[CrossRef]

Thomson, C. J. (2004). Coherent-state analysis of the seismic head-wave problem: an overcomplete representation and its relationship to rays and beams, Geophys. J. Int.157 ,1189 –1205.[CrossRef]

Ulrych, T. J., S. L. M. Freire, and M. D. Sacchi (1999). Eigenimage processing of seismic sections, in Covariance Analysis for Seismic Signal Processing, R. L. Kirlin and W. J. Done (Editors), Society of Exploration Geophysicists, Tulsa, Oklahoma,241 –274.

Vinnik, L. P. (1977). Detection of waves converted from P to PV in the mantle, Phys. Earth Planet. Interiors15 ,39 –45.[CrossRef]

White, B. S., A. Norris, A. Bayliss, and R. Burridge (1987). Some remarks on the Gaussian beam summation method, Geophys. J. R. Astr. Soc.89 ,579 –636.

Wu, R. S., and L. Chen (2001). Beamlet migration using Gabor-Daubechies frame operators, EAGE 63rd Annual Meeting, Amsterdam, p. 74.

Wu, R. L., and L. Chen (2002). Wave propagation and imaging using Gabor-Daubechies beamlets, in Theoretical and Computational Acoustics, E. C. Chang, Q. Li, and T. F. Gao (Editors), World Scientific, New Jersey,666 –670.

Xu, S., and G. Lambere (1998). Maslov+Born migration/inversion in complex media, in 68th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts,1704 –1707.

Yu, J., L. J. Katz, F. Followill, H. Sun, and G. T. Schuster (2003). Autocorrelation migration: IVSPWD test, Geophysics 68,297 –307.[CrossRef][ISI][GeoRef]

Zhang, Y., S. Gray, and J. Young (2000). Exact and approximate weights for Kirchhoff migration, in 70th Ann. Internat. Mtg. Soc. Explor. Geophys., Expanded Abstracts,1036 –1039.

Zhu, L. (2000). Crustal structure across the San Andreas Fault, Southern California from teleseismic converted waves, Earth Planet. Sci. Lett. 179,183 –190.[CrossRef]




This article has been cited by other articles:


Home page
Bulletin of the Seismological Society of AmericaHome page
S. Dasgupta and R. L. Nowack
Deconvolution of Three-Component Teleseismic P Waves Using the Autocorrelation of the P to SV Scattered Waves
Bulletin of the Seismological Society of America, October 1, 2006; 96(5): 1827 - 1835.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF)
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Nowack, R. L.
Right arrow Articles by Sheng, J.-M.
Right arrow Search for Related Content
GeoRef
Right arrow GeoRef Citation


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS