high fidelity seismic imaging and parameter - OSTI.gov

high fidelity seismic imaging and parameter - OSTI.gov

HIGH RESOLUTION/HIGH FIDELITY SEISMIC IMAGING AND PARAMETER ESTIMATION FOR GEOLOGICAL STRUCTURE AND MATERIAL CHARACTERIZATION Final Technical Report ...

8MB Sizes 0 Downloads 7 Views

HIGH RESOLUTION/HIGH FIDELITY SEISMIC IMAGING AND PARAMETER ESTIMATION FOR GEOLOGICAL STRUCTURE AND MATERIAL CHARACTERIZATION

Final Technical Report Period Covered April 15, 2004 – April 14, 2008

Ru-Shan Wu and Xiao-Bi Xie UNIVERSITY OF CALIFORNIA AT SANTA CRUZ Santa Cruz, CA 95064

Submitted December 17, 2008

PREPARED FOR THE UNITED STATES DEPARTMENT OF ENERGY OFFICE OF BASIC ENERGY SCIENCES Work Performed Under Contract No. DE-FG02-04ER15530

Table of Contents  PUBLICATIONS UNDER THIS PROJECT ............................................................................. 4  ABSTRACT ................................................................................................................................ 6  1. INTRODUCTION .................................................................................................................. 7  PART I. RESOLUTION ANALYSIS FOR SEISMIC IMAGING .......................................... 9  2. SPATIAL RESOLUTION OF SEISMIC IMAGING .......................................................... 10  2.1. Background .................................................................................................................... 10  2.2. Formulation of Spatial Resolution Based On Inverse Theory ....................................... 10  2.3. Angular-Spectral Representation of Point Spreading Function (PSF) .......................... 14  2.4. Resolution Analysis for Interface Scattering ................................................................. 18  2.5. Theoretical Resolution of An Acquisition System (Acquisition Resolution) ................ 21  2.6. Numerical Examples of Actual Resolution for Different Imaging Systems.................. 25  2.7. Conclusions .................................................................................................................... 25  3. SEISMIC RESOLUTION AND ILLUMINATION: A WAVE-EQUATION-BASED ANALYSIS ............................................................................................................................... 28  3.1. Background .................................................................................................................... 28  3.2. Formulation .................................................................................................................... 29  3.3. Numerical Examples ...................................................................................................... 32  3.4. Conclusions .................................................................................................................... 34  PART II. AMPLITUDE FIDELITY OF SEISMIC IMAGING: TRUE-REFLECTION IMAGING .................................................................................................................................... 35  4. OPTIMIZATIONS IN TRUE-AMPLITUDE ONE-WAY PROPAGATORS..................... 36  4.1. Implementations of the true-amplitude one-way wave equation ................................... 37  4.2. Optimizations for different terms in the true-amplitude one-way propagators ............. 39  4.3. Numerical results ........................................................................................................... 41  4.4. Conclusions .................................................................................................................... 45  5. SUPERWIDE-ANGLE ONE-WAY MIGRATION: APPLICATION TO 2004 BP 2D BENCHMARK MODEL .......................................................................................................... 47  5.1. Background .................................................................................................................... 47  5.2. Theory ............................................................................................................................ 47  5.3. Experiments on 2004 BP 2D Benchmark Dataset ......................................................... 48  5.4. Discussions and Conclusions ......................................................................................... 50  6. TARGET ORIENTED FULL-WAVE EQUATION BASED ILLUMINATION ANALYSIS ................................................................................................................................................... 56  6.1 Background ..................................................................................................................... 56  6.2 Formulations for Different Illumination Measurements ................................................. 57  6.3. Numerical Examples for Illumination Calculations ...................................................... 58  6.4. Conclusions .................................................................................................................... 67  7. LOCAL-ANGLE DOMAIN ILLUMINATION IN FREQUENCY DOMAIN FOR FULLWAVE PROPAGATORS......................................................................................................... 68  7.1. Wavefield Decomposition for One-Way Propagators ................................................... 69  7.2. Wavefield Decomposition for Full-Wave Propagators.................................................. 70  7.3. Local-Angle Domain Illumination Results .................................................................... 72 

2

7.4. Discussion About The Influence of Multiples on The Illumination Strength ............... 76  7.5. Conclusions .................................................................................................................... 78  Appendix A Gabor-Daubechies Frame (GDF) Beanlet Decomposition .............................. 78  8. FAST ACQUISITION APERTURE CORRECTION BY BEAMLET DECOMPOSITION ................................................................................................................................................... 80  8.1. Beamlet Decomposition of the Wavefield ..................................................................... 81  8.2. LAD Image and Image Amplitude Correction Factor by LWD Bemlet Decomposition ............................................................................................................................................... 82  8.3. Examples ........................................................................................................................ 84  8.4. Discussion ...................................................................................................................... 93  8.5. Conclusions .................................................................................................................... 96  Appendix A: Gabor-Daubechies Frame (GDF) Beanlet Decomposition ............................. 97  Appendix B: Local Exponential Frame (Lef) Beamlet Decomposition ............................... 97  9. A DEPTH MIGRATION METHOD BASED ON THE FULL-WAVE REVERSE-TIME CALCULATION AND LOCAL ONE-WAY PROPAGATION ........................................... 100  9.1. Background .................................................................................................................. 100  9.2. Methodology ................................................................................................................ 100  9.3. Numerical Examples .................................................................................................... 102  9.4. Discussions and Conclusion ........................................................................................ 105  PART III. IMAGING AND INVERSION .............................................................................. 108  10. GENERALIZED DIFFRACTION TOMOGRAPHY IN HETEROGENEOUS MEDIA WITH FINITE DATA APERTURE....................................................................................... 109  10.1. Backpropagation and the Imaging Principle .............................................................. 109  10.2. Imaging Condition in the Local Angle-Domain and the Local Image Matrix .......... 110  10.3. Deconvolution Filtering in Local Wavenumber Domain and the Generalized Diffraction Tomography ..................................................................................................... 111  10.4. Relation with the Classic Diffraction Tomography in A Homogeneous Background ............................................................................................................................................. 115  10.5. Conclusion ................................................................................................................. 116  REFERENCES ....................................................................................................................... 117 

3

PUBLICATIONS UNDER THIS PROJECT Cao, J., and Wu, R.S., 2005, Influence of Propagator and Acquisition Aperture on Image Amplitude, 75th Annual International Meeting, SEG, Expanded abstracts, 1946-1949. Cao, J., and Wu, R.S., 2006, Study of the influence of propagator amplitude correction on image amplitude using beamlet propagator with local WKBJ approximation, 76th Annual International Meeting, SEG Expanded Abstracts, 2499-2503. Cao, J., and R.S. Wu, 2008, Amplitude compensation for one-way wave propagators in inhomogeneous media and its application to seismic imaging, Communications in Computational physics, 3, 203-221. Fehler, M., L. Huang, R.S. Wu, and X.B. Xie, 2005, Seismic image resolution: numerical investigation of role of migration imaging operator, 75th Annual International Meeting, SEG, Expanded Abstracts, 1870-1873. He, Y.F., and Wu, R.S., 2007, One-return boundary element method and salt internal multiples, 77th Annual International Meeting, SEG, Expanded abstracts, 2080-2084. Jia, X.F., and Wu, R.S., 2007, Imaging steep salt flanks by super-wide angle one-way method, 77th Annual International Meeting, SEG, Expanded Abstracts, 2265-2269. Luo, M., R.S. Wu, and X.B. Xie, 2005, True amplitude one-way propagators implemented with localized corrections on beamlets, 75th Annual International Meeting, SEG, Expanded Abstracts, 1966-1969. Mao, J., and Wu, R.S., 2007, Illumination analysis using local exponential beamlets, Expanded abstracts, 77th Annual International Meeting, SEG, Expanded Abstracts, 2235-2239. Wu, R.S., and Jia, X.F., 2006, Accuracy improvement for super-wide angle one-way waves by wavefront reconstruction, 76th Annual International Meeting, SEG, Expanded abstracts, 2976-2980. Wu, R.S., X.B. Xie, and X.Y. Wu, 2006, one-way and one-return approximations (dewolf approximation) for fast elastic wave modeling in complex media, in: Advances in Wave Propagation in Heterogeneous Earth, Elsevier, 265-322. Wu, R.S., X.Y. Wu, and X.B. Xie, 2006, Simulation of high-frequency wave propagation in complex crustal waveguides using generalized screen propagators, in: Advances in Wave Propagation in Heterogeneous Earth, Elsevier, 323-363. Wu, R.S., X.B. Xie, M. Fehler and L. Huang, 2006, Resolution analysis of seismic imaging, 68th Conference & Technical Exhibition, EAGE, Expanded Abstracts, GO48.

4

Wu, R.S., 2007, Generalized diffraction tomography in inhomogeneous background with finite data aperture, 77th Annual International Meeting, SEG, Expanded Abstract, 2728-2732. Wu, R.S., and Mao, J., 2007, Beamlet migration using local harmonic bases, 77th Annual International Meeting, SEG, Expanded abstracts, 2230-2234. Wu, R.S., Wu, B.Y. and Geng, Y. 2008, Seismic wave propagation and imaging using timespace wavelets, Expanded abstracts, SEG 78th Annual Meeting, 2983-2987. Xie, X.B., R.S.,Wu, M. Fehler and L. Huang, 2005, Seismic resolution and illumination: A wave-equation-based analysis, 75th Annual International Meeting, SEG, Expanded Abstracts, 1862-1865. Xie, X.B., S.W. Jin and R.S. Wu, 2006, Wave equation based seismic illumination analysis, Geophysics, 71, S169-177. Xie, X.B. and R.S. Wu, 2006, A depth migration method based on the full-wave reverse-time calculation and local one-way propagation, 76th Annual International Meeting, SEG, Expanded Abstracts, 2333-2337. Xie, X.B. and H. Yang, 2007, A migration velocity updating method based on the shot index common image gather and finite-frequency sensitivity kernel, 77th Annual International Meeting, SEG, Expanded Abstracts, 2767-2771. Xie, X.B., and H. Yang, 2008, A full-wave equation based seismic illumination analysis method, 70th Conference & Technical Exhibition, EAGE, Expanded Abstracts, P284. Xie, X.B., and H. Yang, 2008, The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis, Geophysics, 73, S241-S249. Zhu, X., and R. Wu, 2008, Imaging diffraction points using the local image matrix in prestack migration, 78th Annual International Meeting, SEG, Expanded Abstracts, 2361-2365.

5

ABSTRACT Our proposed work on high resolution/high fidelity seismic imaging focused on three general areas: (1) development of new, more efficient, wave-equation-based propagators and imaging conditions, (2) developments towards amplitude-preserving imaging in the local angle domain, in particular, imaging methods that allow us to estimate the reflection as a function of angle at a layer boundary, and (3) studies of wave inversion for local parameter estimation. In this report we summarize the results and progress we made during the project period. The report is divided into three parts, totaling 10 chapters. The first part is on resolution analysis and its relation to directional illumination analysis. The second part, which is composed of 6 chapters, is on the main theme of our work, the true-reflection imaging. True-reflection imaging is an advanced imaging technology which aims at keeping the image amplitude proportional to the reflection strength of the local reflectors or to obtain the reflection coefficient as function of reflection-angle. There are many factors which may influence the image amplitude, such as geometrical spreading, transmission loss, path absorption, acquisition aperture effect, etc. However, we can group these into two categories: one is the propagator effect (geometric spreading, path losses); the other is the acquisition-aperture effect. We have made significant progress in both categories. We studied the effects of different terms in the true-amplitude oneway propagators, especially the terms including lateral velocity variation of the medium. We also demonstrate the improvements by optimizing the expansion coefficients in different terms. Our research also includes directional illumination analysis for both the one-way propagators and full-wave propagators. We developed the fast acquisition-aperture correction method in the local angle-domain, which is an important element in the true-reflection imaging. Other developments include the super-wide angle one-way propagator and special full-wave reverse-time migration method. Finally, we studied the theoretical basis of true-reflection imaging and bridges imaging and inversion with the theory of diffraction tomography.

6

1. INTRODUCTION Reflection seismic migration imaging is used to determine the locations and geometric characteristics of structures in the Earth that are of interest in resource identification, environmental remediation efforts, and for basic science applications. Due to the success in obtaining improved images of the geometry of subsurface structures with wave-equation migration imaging, there is now interest in obtaining more information from seismic data than the geometry of subsurface strata. To advance the applicability of wave equation imaging methods for estimating rock and pore fluid properties, our proposed work aimed to develop improved high resolution/high fidelity methods for wave-equation based imaging and to extend these methods for estimating impedance contrasts at boundaries. Our proposed work focuses on three general areas: (1) development of new, more efficient, wave-equation-based propagators and imaging conditions, (2) developments towards amplitude-preserving imaging in the angle domain, in particular, imaging methods that allow us to estimate the reflection as a function of angle at a layer boundary, and (3) studies of wave inversion for local parameter estimation. This report is divided into three parts. Part I is “RESOLUTION ANALYSIS FOR SEISMIC IMAGING”. In this part, first (Chapter 2) we summarize our theoretical study on resolution analysis. We start from the general formulation of resolution operator based on the inverse theory. We applied the local angle-domain decomposition to the resolution operator and derived the angular-spectral representation of the resolving kernel. Our resolution analysis including the effects of both the acquisition system and the imaging process and put the resolution analysis on a sound theoretical basis. Then (Chapter 3) we relate the resolution analysis to the directional illumination analysis. Through theoretical analysis and numerical examples, we demonstrate how the acquisition-system defines the resolution limit and how to improve the resolution by using accurate imaging methods. In part II we discuss “AMPLITUDE FIDELITY OF SEISMIC IMAGING: TRUEREFLECTION IMAGING”. True-reflection imaging is an advanced imaging technology which aims at keeping the image amplitude proportional to the reflection strength of the local interface. True-reflection image can present the image as the total strength or the reflections as a function of reflection-angle. In our study, we mainly concentrated on the total strength imaging. There are many factors which may influence the image amplitude, such as geometrical spreading, transmission loss, path absorption, acquisition aperture effect, etc. However, we can group these into two categories: one is the propagator effect; the other is acquisition-aperture effect. In Chapter 4 we discuss the wave-equation based true-amplitude propagator. We discussed the different terms in the true-amplitude one-way propagator, especially the terms include lateral velocity variation of the medium. We also demonstrate the improvements by optimizing the expansion coefficients in different terms. In Chapter 5, we introduced a new super-wide angle one-way propagator, which expands the propagate angle to beyond 90 degrees and therefore can model turning waves and reflection waves. The application example for the BP model showed its excellent property. In order to overcome the limitation of one-way wave propagator, there is a trend in the industry to use the full-wave propagator for seismic imaging, even with significantly increased cost. To apply the concepts of true-reflection imaging to the full-wave imaging methods, such as

7

the finite difference reverse-time migration, we need to analyze the acquisition-aperture effects using directional illumination analysis. In Chapters 6 and 7 we summarize our results on directional illumination analyses using time-space domain method (Chapter 6) and frequencydomain method. These methods can serve as illumination analysis during acquisition design, quality control of the imaging system, and the image amplitude analysis and compensation of full-wave migration such as reverse-time migration. In Chapter 8 we present the results of our study on fast acquisition-aperture correction method in true-reflection imaging. The method uses the beamlet decomposition developed by our group to the wavefield and Green’s function decompositions into the local wavenumber domain, resulting in significant increase of computation efficiency (more than twice faster) for acquisition-aperture correction. This paved the way for practical application of true-reflection imaging in the exploration industry. In Chapter 9 a depth migration method of full-wave reverse-time migration which uses local one-way propagation to implement the imaging condition is presented. In Part III “IMAGING AND INVERSION” we present the preliminary results of our study on diffraction tomography in heterogeneous media with finite data aperture. This study bridges the gap between imaging and inversion. It takes the imaging process as the first step of a local inversion. The theory is based on the resolution operator and generalizes the classic diffraction tomography to the case of heterogeneous media with finite data aperture (including the spatial aperture and the frequency band). It also takes the true-reflection imaging as a special case of the tomographic inversion.

8

PART I. RESOLUTION ANALYSIS FOR SEISMIC IMAGING

9

2. SPATIAL RESOLUTION OF SEISMIC IMAGING 2.1. Background Spatial resolution has been studied by Beylkin et al. (1986) based on generalized Radon transform and the mapping of domain integration (frequency band and spatial aperture) into a spectral coverage in spatial frequency domain (wavenumber domain). The mapping is done by simple ray-tracing which does not take into account of the frequency dependence and other wave phenomena. Since then, the topic has been investigated by many authors from different points of view (Chen and Schuster, 1999, 2001; Gelius et al., 2002; Gibson and Tzimeas, 2002). While these analyses are useful tools for the resolution problem, most analyses are formulated to calculate the influences of acquisition system to resolution, and did not look at other factors in the imaging/migration process, such as the accuracy of the propagators, errors in velocity model, etc. Some authors called such a resolution as the resolution limit, the best resolution of an acquisition system can get. In this research, we formulate the resolution of imaging system, which includes the acquisition system and imaging (migration) process. The analysis is based on the inverse theory and the local angle domain decomposition of Green’s functions. This resolution of imaging system represents and quantifies the actual resolution we observed in the migrated images and provides a theoretical basis for the estimation of different factors influencing the resolution and image quality in complex media. The other important feature of our formulation is that the resolution problem is treated totally based on wave theory and no high-frequency asymptotic approximation is made.

2.2. Formulation of Spatial Resolution Based On Inverse Theory Spatial resolution can be studied under the general frame of inversion theory. Resolution operator or the discrete form, resolution matrix has been defined to quantify the resolution of parameter inversion by a particular inversion scheme (Bucks and Gilbert, 1070; Aki and Richards, 1980, section 12.3; Tarantola, 1987, section 4.33, 7.22). The resolution matrix has dependence on both the acquisition system and the inversion scheme. Assume the acquisition process can be modeled by Fm  d (2.1) Where F is the forward modeling operator and d is the data. If we adopt a specific inversion operator B to apply to the data we can get a set of model parameters mI, which is different from the real m, mI  Bd (2.2) Substituting (2.2) into (2.1) we get the relation between the inverted model and the real model mI  BFm (2.3) and the resolution matrix (or operator) is defined as R  BF (2.4) For an exact inversion of a well-posed problem, we should have RI (2.5)

10

where I is the identity matrix. In general, the resolution matrix is not an identity matrix and the spreading of the matrix elements along the diagonal and off-diagonal give some quantitative measure of the parameter resolution of the inversion. For the sake of simplicity, here we will use the noiseless formulation. For the stochastic approach (Tarantola, 1987) a similar derivation can be obtain. The imaging problem can be formulated as a specific inverse problem. We assume the media can be decomposed into a smooth variation of velocity and sharp jumps of impedance (discontinuities). The smooth part of heterogeneities can be treated as background (reference) media and the sharp boundaries are left for the inversion (imaging). The velocity distribution of the background media can be derived with different approaches and is assumed known in the “imaging problem”. The unknowns in the imaging problem are the parameter strengths and their locations (distribution). For the problem of spatial resolution, we assume the scattering coefficients are known everywhere. Therefore the resolution matrix is totally defined by the spatial resolution. In the following, we will use the migration operator (simply backpropagation integral) plus the zero-time imaging condition (Claerbout, 1985) as the inversion operator (here the imaging operator) to show the effects of different factors to the final resolution. We will also discuss the influence of amplitude corrections (true-amplitude imaging) to the resolution in later sections. First we discuss the case of volume scattering for scalar waves, where the scattering potential s ( x) is a scalar field with an isotropic scattering pattern. We can write the space-domain formulations for modeling (acquisition process) and imaging (migration) as:  xs  As  (2.6) us ( , xg ; xs )   k 2  dx ' GM ( , x '; xs ) s ( x)GM ( , xg ; x ')  xg  Ag ( xs ) V   A f 

I (x)   d  dxsWs ( , xs )GI ( , x; xs )2  dxgWg ( , xg ) *

Af

As

Ag

GI * ( , x; xg ) z

us ( , xg , xs ) (2.7)

where us is the scattered wave field observed on the surface at xg excited by a source on the surface at xs . GM is the Green’s function of the modeling (acquisition) process, which may include all the factors (geometric spreading, intrinsic and scattering attenuation, boundary scattering, etc.) for the real heterogeneous media; GI is the Green’s function of inversion process, which could be quite different from GM . The integrations on the receiver aperture and source aperture are for the prestack migration process, W’s are the weighting functions for integrations. The integration on the receiver aperture Ag is the Rayleigh integral which simulates the backpropagation process. The integration on frequency is from the imaging condition which states that the downward extrapolated source field and the scattered field will meet at zero time at the scattering points. In (2.7) we used the cross-correlation imaging condition. Other imaging conditions for correcting the imaging amplitude can be also applied, but the general conclusion of resolution analysis should remain unchanged.

11

Equation (2.6) is formulated for volume perturbations (heterogeneities), such as in the case of tomographic inversion or other parameter inversion. In the formulation, s ( x) is assumed as a scalar quantity, and the scattering pattern due to this point scatterer is isotropic. However, in many cases the local scattering pattern is not isotropic, such in the case of elastic wave scattering, or boundary scattering. Especially in seismic migration/imaging, it has been accepted in the community that the migration/imaging process is a discontinuity (or boundary) imaging (Beylkin, 1985; Bleistein et al, 2001). If we can approximate the curved boundaries with assemble of linear boundary elements, then for each boundary element the local scattering potential becomes angledependent s (x, n ) , where  n is the dip-angle of the boundary element. More generally, the scattering potential may become s (x, xs , xg ) . In this way the local scattering property can be represented by a local scattering matrix (LSM) (Wu and Chen, 2002; Wu et al, 2004), in which both the incident and scattering angles are independent variables. This angle dependent scattering can be also formulated from the theory of scattering, and will be discussed later. For now let us consider the case of scalar potential (point scatterer with isotropic scattering pattern). Write modeling and imaging process (2.6) and (2.7) into operator form, resulting in U( , xg , xs )  F ( , xg , xs | x 0 )S(x 0 )

(2.8)

I (x)  B(x |  , xs , xg )U( , xs , xg )

(2.9)

where F is the acquisition (modeling) operator which maps the model S into the data set U; while B is the imaging operator which invert the data U into the image I. Therefore the resolution operator is obtained as R (x; x 0 )  B(x |  , xs , xg )F ( , xg , xs | x0 )

(2.10)

The kernels for the resolution operator can be obtained as (x; x0 )    d k 2  dxsWs ( , x, xs )GI * ( , x; xs )GM ( , x 0 ; xs ) As

2  dxgWg ( , x, xg ) Ag

GI * ( , x; xg ) z

12

GM ( , xg ; x 0 )

(2.11)

Figure 2.1. Schematic diagrams showing the modeling (data acquisition) and imaging (inversion) processes. For zero-offset (zero receiver aperture) acquisition, or monostatic measurement in radar terminology, the above defined resolution matrix is degenerated to GI * ( , x; xs ) 2 * (x; x0 )  2 dk  dxsWs ( , x, xs )GI ( , x; xs )GM ( , x0 ; xs ) GM ( , xs ; x 0 ) As z (2.12) Note that the amplitude functions in (2.12) are different from the exploding-reflector modeling. In the latter case, (2.12) is further simplified to (x; x 0 )  2  d k 2  dxgWg ( , x, xg ) Ag

GI * ( , x; xg ) z

GM ( , xg ; x 0 )

(2.13)

The above derived resolution matrix (operator) is in fact the spatial resolution matrix for the whole acquisition and migration process (the imaging system). The matrix element (x0 , x) is called the resolving kernel of the resolution operator (Backus and Gilbert, 1970: Tarantola, 1987), which is the point spreading function of the imaging system. In this way the point spreading function (impulse response) is defined under the guidance of general inversion theory. If we set the scatterer’s distribution s (x ')   (x ' x0 ) in equation (2.6) and substitute it into equation (2.7), we see the equivalence of the imaging process to the calculation of resolution matrix. This gives us the numerical procedure of calculating the resolution matrix or PSF for any imaging system (including acquisition and migration)/inversion). In the same way the PSF (point spreading function) of zero-offset imaging process or imaging with exploding-reflector modeling defined in (2.12) and (2.13) respectively, can be calculated by numerical simulations. We see that the point spreading function (PSF) is dependent on both the data acquisition and the migration (imaging) process (including the propagators and the imaging condition). The Green’s function GM in the modeling integral is the real Green’s function in the scattering process, which may include absorption and scattering losses and geometric spreading. However,

13

the GI in the migration (backpropagation) process can be quite different depending on the backpropagation operator. In fact, there are so many different propagators, e.g. full-wave FD propagators, one-way FD propagators, one-way dual-domain propagators, ray propagators, are used in the migration practice, and they will have different resolution matrices and PSF’s. The PSF from the acquisition system itself can only provide the theoretical resolution, or the resolution limit.

2.3. Angular-Spectral Representation of Point Spreading Function (PSF) Angular-spectral representation of resolution of PSF is more intuitive. We can see directly the information coverage in local angle-domain. We perform local 3D Fourier transform on (x, x0 ) with coordinate center at x 0 to get its 3D spectrum: (K ; x 0 )   dxe  iK  ( x  x0 )W (x, x 0 )(x; x 0 ) v

(2.14)

where K is the 3D local wavenumber vector and K  K  , with  as the unit angular vector, W (x, x0 ) is a 3D window function to localize (x; x0 ) . For simplicity, we set the weighting functions Ws and Wg in migration as unity. It is known that the window spectrum W (K , x0 ) will convolve with the obtained PSF spectrum. Therefore the window W (x, x0 ) acts as a smear operator in spatial frequency domain, resulting in the reduction of angular resolution of the PSF spectrum. For resolution analysis based on wave theory, the spectral resolution and the spatial localization have to satisfy the uncertainty principle. It means that if you want more localized PSF (resolving kernel) in space domain to study the PSF change in strongly heterogeneous media, you have to sacrifice a little more on directional localization (spectral resolution). Substituting (2.11) into (2.14) we obtain (K ; x 0 )    d k 2  d K s  dxs GI * ( , x 0 , K s ; xs )GM ( , x 0 ; xs ) Af

2  dxg Ag

As

GI * ( , x 0 , K  K s ; xg ) z

(2.15) GM ( , xg ; x 0 )

where G ( , x 0 , K ; x ') is the beamlet decomposition (or local plane-wave decomposition) of the Green’s function around x0 , and x’ is the source or receiver position. The 3D wavenumber K s is related to the propagating wavenumber through K s  ( s ,  s )  k s ,

(2.16)

Where  s is the horizontal wavenumber and

 s  k 2   s2

14

(2.17)

is the vertical wavenumber with k   / v( x 0 ) , and v( x 0 ) as the local velocity. In (2.15) Af is the frequency aperture (band), and the integration over K s is a convolution integral in K domain. Therefore, K g  K  K s  ( g ,  g )  k g (2.18) The spectral component determined by this scattering experiment with a single frequency will be at K  kg  ks  k (2.19) where k  is the exchange wavenumber. Therefore, (2.15) can rewritten as (K  k  ; x 0 )    d S ( )k 2  dxs GI * ( , x 0 , k s ; xs )GM ( , x 0 ; xs ) Af

2  dxg Ag

As

GI * ( , x 0 , k g ; xg ) z

(2.20)

GM ( , xg ; x 0 )

where S ( ) is the source spectrum. If spectral whitening has been incorporated into the imaging process, then we have S ( )  1 . Beamlet decomposition decomposes the field of a Green’s function, a spherical wave front, into beamlets. Each beamlet has a center location x0 with certain width, and a lobe in the specified direction. Because the existence of the lobe, the directional resolution is not as sharp as in the global Fourier transform. The spatial localization of the beamlet is not a point, but with certain width. The compromise between the spatial and directional localizations is governed by the uncertainty principle. In (2.20), GI * ( , x 0 , k s ; xs )GM ( , x 0 ; xs ) can be considered as a point source at x 0 radiating wavefield to x s , and then backpropagated (phase conjugate) again to x 0 , with the local plane wave k s . If we use the high-frequency asymptotic approximation of Green’s function, then the incident wave direction will be determined by the ray direction. For the more general case of wave propagation, the local plane wave decomposition is more appropriate (Wu and Chen, 2002; GI * ( , x 0 , k g ; xs ) GM ( , xg ; x 0 ) Xie and Wu, 2002; Wu et al. 2004). In a similar way z corresponds to a point source at x 0 , radiating wavefield to x g , and then backpropagated to x 0 , with the local plane wave k g . The exchange wavenumber and the related angles can be calculated as

k   2k sin( sc / 2) , and

 K  ( g   s ) / 2,

 r   sc / 2  ( g   s ) / 2 ,

(2.21)

(2.22)

where  K is the unit direction vector of K  k  , k   / c0 , with c0  v( x 0 ) and  sc is the scattering angle defined as the angle between the incident direction and the receiving direction.

15

 K is the angle of K vector with respect to a reference direction, such as the z-axis;  s and  g are the angles of the k s and k g vectors respectively, and  r is the reflection angle when dealing with local reflectors. For each frequency, we may find a pair of (k s , k g ) satisfying K  k s  k g . The final magnitude of the spectral component (K ; x0 ) is the sum of contributions from all the pairs. The effect of data aperture (domain of integration), including the frequency band and acquisition spatial aperture, is in the integration apertures Af , As and Ag . Equation (2.20) also provides a numerical implementation for the calculation of angular spectra of resolution matrix. We need first to decompose the Green’s function GI * ( , x; xs ) , GI * ( , x; xg ) / z into the local angle-domain (local source and receiver angles) at x0 :

GI * ( , x0 , k s ; xs ) , GI * ( , x 0 , k g ; xg ) / z . The decomposition can be done using beamlet

transform (window Fourier frame) (Wu and Chen, 2002; Wu et al. 2004) or local slant stack (Xie and Wu, 2002). In this way, we can calculate the resolution (PSF) for any points in a complex model. If the angle-domain Green’s functions can be stored for the whole model, we can have an efficient algorithm for analyzing the resolution for any target point in the image space.

Figure 2.2: Spectral representation of the resolving kernel (point spread function) for volume scattering at a point in a normal area (point A) of the SEG/EAGE salt model: on the left is the amplitude spectrum, and on the right, the phase spectrum

16

Figure 2.3: Same as Figure 2.2 but for a point in the subsalt area (point B) Figure 2.2 and 2.3 give the spectral representation of the resolution at two selected points in the SEG/EAGE salt model. One point is in a normal region (point A in Figure 2.2) and the other point in the subsalt region (point B in Figure 2.3). We see that due to the influence of the irregular salt-body and the acquisition configuration, the angular spectrum of PSF is quite nonuniform and the resolution lengths are very different for different directions. In the similar way we can write the PSF in angular spectrum domain for zero-offset and exploding reflectors modeling respectively: (K  k  ; x 0 )  2 dk 2  dxs GI * ( , x 0 , k s ; xs )GM ( , x 0 ; xs ) As

(K  k s ; x0 )  2 dk 2  dxs As

GI * ( , x 0 , k g ; xg ) z

GM ( , xs ; x0 ) , (2.23)

GI * ( , x 0 , k s ; xs ) GM ( , xs ; x0 ) . z

(2.24)

Under asymptotic approximation, we can put k s  k g and K  k   2k s in (2.23). Directional resolutions:

It is convenient to study the directional resolution using the spectral representation of the resolving kernel (2.20). Due to the correspondence between the spatial and spectral distributions, the resolution length along certain spatial direction is inversely proportional to spectral width along that direction. Based on this, we can define the maximum resolution, depth and horizontal resolutions, vertical and lateral resolutions, etc. 17

2.4. Resolution Analysis for Interface Scattering Many kinds of heterogeneous media in the earth can be decomposed into a rough part and a smooth part. The rough part is composed of irregular sharp interfaces. The smooth part is formed by many smoothly varying domains surrounding by those sharp interfaces. The smoothly varying parts of the medium cause diffraction, path bending, travel-time variation, amplitude changes of the wavefield (forward scattering), but only sharp boundaries produce significant backscattering (reflection) and large-angle scattering. In this case, we can formulate the modeling process with boundary scattering model, in which the scattering only occurs at the sharp interfaces, so that the volume integral for scattered waves in (2.6) becomes a boundary integral. The boundary scattering can be calculated with the Kirchhoff integral: us (x)   us (x ') 

  G (x, x ')ds (x ')   G (x, x ') u (x ')ds (x ') ,  n n s

(2.25)

where ds (x ') is a boundary element at x ' on the boundary  , and n is the outward normal at that point. In order to calculate the scattered field and its normal gradient on the interface, we apply the Kirchhoff approximation, which is also called the physical optics approximation or tangent plane approximation (Ogilvy, 1991). At any point of the curved interface, based on the tangent plane approximation we can define the local reflection coefficient as R(x ', nˆ , r ) , where x ' is the location on the surface, nˆ is the outward normal unit vector at the point, and  r is the reflection angle with respect to the normal. In the case of surface reflection profiling on the heterogeneous earth, various one-way propagators can be used as the Green’s function. The source and receiver (geophone) fields can be downward propagated to a level close to the curved interface (Figure 2.4), and then decomposed into local plane waves (Wu and Chen, 2002, 2006; Xie and Wu, 2002). In fact the local plane waves derived by local windowing are beamlets with a dominant direction and certain lobes (ibid). In order to define the local reflection coefficients, we need also to decompose the Greens functions into local plane waves at depth z. Then equation (2.25) becomes u s ( x g , x s )  i   s g (2.26)  GM ( s , x ', z; xs )GM ( g , x ', z; x g ) R(x ', n垐,r )[(k g  k s )  n]exp{i(k g  k s )  r '}ds(x ') 

where GM ( s , x ', z; x s ) is the coefficients of a local plane wave component of the Green’s

function, b jl  x  is beamlet atom (element function), and g  x  xl  is the window function for decomposition (spatial localization), and  j is the local wavenumber (see Wu and Chen, 2002,

2006). In above equation, r '  ( x ', z ' z ) , and k s  ( s ,  s ) is the local wavenumber vector, with

 s  ( / v) 2   s 2 as the vertical wavenumber, for the local plane wave component. It is seen that the original space-domain Green’s function is transformed into a beamlet domain Green’s function (Green’s beamlets), which represents the local plane wave response near the boundary generated by a point source on the acquisition surface. The same decomposition can be applied

18

to the receiver’s Green’s function and GM ( g , x ', z; x g ) is its local plane wave coefficient. For the detailed derivation of (2.26) see Wu (2006).

Figure 2.4: Local plane wave decomposition near the curved interface Similar to (2.20) for the case of volume scattering, we obtain the resolving kernel of spatial resolution (PSF) in the local wavenumber domain for interface scattering: (K  k  ; n垐 , x ',  r )   d S ( )  dxs [k   n]GI* (x ', k s ; xs )GM (x ', k s ; xs ) Af

2  dxg Ag

As

GI* (x ', k g ; xg ) GM (x ', k g ; xg ) z

.

(2.27)

where S ( ) is the source spectrum. Figure 2.5 shows the definition of the local wavenumbers and the asymptotic sphere for reflector imaging by a single-frequency. For a given reflection angle  r (the corresponding magnitude of k  is 2k sin  r ), the resolving kernel defines how the imaging system maps a single dip into a distribution on the asymptotic sphere with k   k g  k s , in the local angle-domain at the image point. Comparing (2.27) with the resolving kernel for volume scattering (2.20), we can see the difference between the two cases. The spectral coverage for boundary scattering is defined on a h-f asymptotic sphere (a circle in the 2D case). The resolution is only for the dip-direction (lateral resolution) and do not have resolution in the magnitude of the spectral component. This is due to the assumption of parameter discontinuity at the boundary reflection point. For a single-dip boundary element with normal nˆ , the spreading pattern in the asymptotic circle determines the lateral resolution of the boundary element. Figure 2.6 shows the PSF for different boundary reflection elements R(x ', nˆ , r ) at point A outside the salt body. Figure 2.7 gives that for point B in the subsalt area. In the figures, the horizontal coordinate is angle of k  of the imaging system, and along the vertical coordinate we marked the real dips for each PSF curves. We can see the different angular spreading for the given R(x ', nˆ , r ) . The range of detectable dip-angles can be studied by the directional illumination analysis. The PSF in spectral domain is closely related to the spectral illumination of the acquisition system. We will discuss this relation further in the next section. For multi19

frequency imaging, the superposition of asymptotic spheres at different frequencies will define the vertical resolution, which essential is the source wavelet, or a pulse proportional to the source bandwidth in the case of spectral whitening.

Figure 2.5: the definition of the local wavenumbers and the asymptotic sphere of reflector imaging by a single-frequency

Figure 2.6: PSF for a boundary element at point A (outside the salt body of the SEG/EAGE salt model) with different dips: (a) reflection-angle is fixed at 0o; (b) reflection angle is fixed at 20o.

Figure 2.7: PSF for a boundary element at point B (in the subsalt area of the SEG/EAGE salt model) with different dips: (a) reflection-angle is fixed at 0o; (b) reflection angle is fixed at 20o. 20

2.5. Theoretical Resolution of An Acquisition System (Acquisition Resolution) If we assume that an exact Green’s function GM is used as the inverse propagator GI , then the effect of imperfect propagator can be eliminated, and the resolution is totally determined by the data aperture. From (2.20) we see that around x0 , the phase of GI will be totally compensated by G *M . Therefore it is sufficient to investigate only the amplitude spectra in the resolution analysis. max (K  k  ; x 0 )   d k 2  dxs GI ( , x 0 , k s ; xs ) GM ( , x 0 ; xs ) Af

2 dxg Ag

As

GI ( , x0 , k g ; xg ) z

(2.28)

GM ( , xg ; x 0 )

The resolution derived this way is a theoretical limit of the acquisition system, i.e. the maximum resolution the imaging system can achieve with the give acquisition configuration. In the above formula, we see that except for a geometric spreading factor GM , it is equivalent to the formulation of resolution analysis in the literature (e.g. Beylkin et al., 1986; Gelius et al., 2002). The calculation of (2.28) can be simplified by some approximations based on the energy conservation principle. The back propagation integral in (2.20) represents a refocusing process of the forward propagated beamlet (spreading) wavefront. For large enough receiver (or source) aperture, we consider that all the energy collected by the receiver array can be refocused to the lobe of the beamlet. Then energy conservation exists GI* ( , x 0 ; xg ) (2.29) |  dxg GM ( , x 0 , k g ; xg ) |2   dxg | GM ( , x 0 , k g ; xg ) |2 Ag Ag z Same reasoning can be applied to the source array back propagation integral if we replace the Green’s function GI ( , x0 ; xs ) by its vertical gradient GI ( , x0 ; xs ) / z in (2.20). Then, we can simplify the calculation of (2.28) to max (K  k  ; x0 )   dk 2   dxs GI ( , x0 , k s ; xs )  As Af

2

1/ 2

2  Ag dxg GI ( , x0 , k g ; xg ) 

.

(2.30)

which is similar to the calculation of the acquisition dip-response (ADR) in the directional illumination analysis (Wu and Chen, 2002, 2006; Xie and Wu, 2002; Xie et al., 2006). However, we need to keep in mind the approximation we made to reach the simplification. Similarly, we can obtain the acquisition resolution for interface scattering model based on (2.27), , x ',  r )   d S ( )   dxs [k   n] GI (x ', k s ; xs ) max (K  k  ; n垐  As Af

21

2

2 Ag dxg GI (x ', k g ; xg ) 

1/ 2

(2.31)

As we stated above, the actual resolution should be also dependent on the accuracy of the backpropagation and reconstruction. This point will be demonstrated further later using propagators with different approximations.

Figure 2.8: Comparison of the “imaging” resolution (the left panel) using the GSP one-way wave imaging propagator to the “acquisition” resolution (the right panel), which gives the resolution limit of the acquisition system: On the top is the resolution at the normal point (point A) of the SEG/EAGE salt model in Figure 2.2; on the bottom is that at the subsalt point (point B) in Figure 2.3.

Figure 2.8 shows the comparison between the “imaging” resolution (the left panel) calculated by equation (2.20) and the “acquisition” resolution (the right panel) calculated by equation (2.30). For the imaging resolution, the modeling Green’s function is computed with a finite difference algorithm, and the imaging (migration) Green’s function is a GSP (generalized screen propagator) one-way wave propagator. For the acquisition resolution, only the one-way propagator is used. On the top of the figure is the resolution at the normal point (point A) of the SEG/EAGE salt model in Figure 2.2; on the bottom is that at the subsalt point (point B) in Figure 2.3. We see that although the basic features are the same, there are some differences between these two resolutions, especially if we consider the phase spectra in Figure 2.2 and 2.3. The similarity between the imaging resolution and acquisition resolution in this case, is due to the exact velocity model used and the good accuracy of the one-way propagator.

22

Figure 2.9 shows the single-frequency acquisition resolution of the acquisition system for the case of SEG/EAGE salt model. The propagator is the one-way dual domain propagator (generalized screen propagator). Although the amplitude may be different from the case of using the actual propagator in data acquisition, the general shapes of the spectral coverage should be close to the theoretical limit of the acquisition system. We see that for the point in the smooth background (the upper left point) the resolution spectral coverage is quite broad. On the other hand, the spectral coverage of the two points in the subsalt region is distorted in different degrees due to the influence of the salt body to the angular coverage. Figure 2.10 gives the acquisition resolution of multi-frequency acquisition. Due to the broad frequency band of the source spectrum, the spectral coverage is significantly improved, especially for low wavenumber components. In comparison, Figure 2.11 gives the theoretical PSF (acquisition resolution) of multi-frequency acquisition (f = 6-24 Hz) for a reflection element with dip-angle 45 degrees for a common scattering-angle gather imaging (  r  20o ). To better understand the features of these angular spectra of resolving kernels, we need to take into consideration of the acquisition configuration for generating the synthetic dataset of the 2D SEG/EAGE salt model. The data acquisition used in the data generation is a simulation of the off-shore acquisition using streamer receiver lines with shot in the right-hand side of the streamer.

Figure 2.9: Single-frequency theoretical PSF (point spreading function) in the angular spectrum domain at some points for the SEG/EAGE salt data acquisition system.

23

Figure 2.10: Multi-frequency theoretical PSF (point spreading function) in the angular spectrum domain at some points for the SEG/EAGE salt data acquisition system.

Figure 2.11: Theoretical PSF (acquisition resolution) of multi-frequency acquisition (f = 6-24 Hz) for a reflection element with dip-angle 45 degrees for a common scattering-angle gather imaging (  r  20o ) at two points of the SEG/EAGE salt model. 24

Definition of resolution lengths based on the PSF: Depth resolution and horizontal resolution: From the PSF (resolution matrix in spatial domain), we measure the resolution lengths in depth direction and in horizontal direction (half value points). Alternatively, we can measure the corresponding widths in wavenumber domain. The widths of spectral coverage in these directions are inversely proportional to the spatial resolution lengths. Vertical resolution and lateral resolution for boundary reflection/scattering elements: For boundary scattering elements, we can measure the corresponding widths along the interfacenormal and transversal directions. In the same way, we can measure the spectral coverage along the normal and transversal (lateral) directions. These resolution lengths are more representative for boundaries with certain dips.

2.6. Numerical Examples of Actual Resolution for Different Imaging Systems The actual resolution will include the effects from the imaging (migration) process expressed by the backpropagation and focusing (imaging condition) in addition to the acquisition effects. To show the difference between the real resolution for an imaging system and the theoretical resolution, we give an example of imaging using exploding reflector data in random media. The random media is given in Figure 2.12a with the exploding point marked by a star. Both one-way migration and ray-Kirchhoff migration are used in the imaging process and the images are shown in Figure 2.12b, c and d, respectively. Since the medium has random heterogeneities with scale comparable to the wavelength of the central frequency, the ray approximation of the Green’s function can produce large errors for long range propagation. Figure 2.13 shows the comparison between the angular spectra of the theoretical resolution (acquisition resolution) (on the top) and the actual resolution in the migrated image (image resolution) (on the bottom). We see the distortion and shrink of the spectrum due to the inaccurate phase information of the raypropagator (bottom right). In comparison, the spectrum given by the one-way wave-equation propagator has much less distortion. Figure 2.14 shows the comparison of lateral resolutions of different propagators.

2.7. Conclusions The resolution of an imaging system defined in this paper based on the inverse theory and the local angle domain decomposition of Green’s functions, includes both the effects of acquisition system and imaging (migration) process. It represents and quantifies the actual resolution we observed in the migrated images. If we neglect the factors such as the errors in backpropagation (migration), and consider only the influence of the acquisition system (frequency band and spatial aperture), the resolution obtained is the resolution limit of the system, the best resolution an acquisition system can get. Theoretical analysis and numerical examples have shown the importance of propagator accuracy to the imaging resolution in complex media.

25

Figure 2.12: Point spreading functions (PDF) of wave equation imaging and ray-Kirchhoff imaging in a random media with small scale heterogeneities.

Figure 2.13: Comparison of the theoretical resolution (acquisition resolution) and the actual resolution in the spectral domain (wavenumber domain): On the top is the angular spectrum of the PSF (point spreading function) calculated from acquisition resolution; on the bottom left is

26

the actual PSF (amplitude spectrum) calculated from the wave-equation imaging; on the bottom right is the PDF from the ray-Kirchhoff imaging.

Figure 2.14: Lateral resolutions at depth 2 km (see Figure 2.7) of wave-equation imaging and ray-Kirchhoff imaging in a random medium.

27

3. SEISMIC RESOLUTION AND ILLUMINATION: A WAVEEQUATION-BASED ANALYSIS 3.1. Background The effect of acquisition on resolution of a seismic image can be evaluated using the coverage of scattering wavenumbers in the target Fourier space (Beylkin, et al., 1985). This coverage is affected by the acquisition configuration, the background velocity model and the frequency band of the signal. Due to the irregular acquisition geometry and complex velocity model involved in the seismic migration imaging processes, the resolution analysis is usually conducted using simplified model geometries or using the ray-based high-frequency asymptotic methods (e.g., Schuster and Hu, 2000; Gibson and Tzimeas, 2002; Yu and Schuster, 2003). The seismic illumination analysis shares many common basis with the resolution analysis (Muerdter and Ratcliff, 2001ab; Muerdter et al., 2001; Berkhout et al., 2001; and Volker et al., 2001). However, traditionally illumination analysis mostly focuses on how the model space is covered by seismic energy and there is no emphasis on image distortion caused by illumination. Recently, progresses have been made in several related areas. The angle related information is emphasized in the illumination analysis and the relationship between the illumination and resolution being investigated (Gelius, et al., 2002; Lecomte, et al. 2003). The local angle related information has been extracted from the wave-equation-based propagators and applied to the illumination analysis (Xie and Wu, 2002; Wu and Chen, 2002, 2003; Wu et al. 2003; Xie et al. 2003, 2004). In this paper, we investigate the relationship between the resolution and illumination using the wave-equation-based propagator. The resolution function can be calculated in both wavenumber and space domains. The new method has the advantage of handling wave motion phenomena. The complex models can be treated without smoothing velocity. Finite frequency band can be used in the calculation. This method can also handle irregular acquisition geometries. To demonstrate potential applications of this method, numerical examples are calculated using the 2D SEG/EAGE salt model. This analysis neglects the influence of different approximations of migration operators and variations of the model from the true structure to the resolution. For a full theory of imaging resolution including the influence of imperfect propagators, and the related numerical examples, see our other studies (Wu, et al., 2005; Fehler, et al., 2005).

Figure 3.1. Coordinate system used in the analysis.

28

3.2. Formulation Consider using a survey system composed of a source located at

rs

and a receiver located at

rg

to

investigate the subsurface target within a small region V  r  neighboring location r (see Figure 3.1). With one-way Green’s functions for heterogeneous background velocity model, the reflection seismic data can be expressed as u (r, rs , rg )  2k02 V m  r   G  r ; rs G  r ; rg  dv , (3.1) where m  r    c c  r   is the velocity perturbation, k0   c0  r  is the background wavenumber and c0 is a local background velocity. Equation (3.1) and the following equations are in frequency domain but we omit the apparent frequency variable. The prestack depth image I at location r  within V can be expressed as



 

  2k02 G*  r ; rs  G*  r ; rg  V m  r   G  r ; rs G  r ; rg  dv .

I (r, r , rs , rg )  G*  r ; rs  G* r ; rg u r, rs , rg

The Green’s functions G in above equations can be decomposed into plane waves within









ik r 

G r ; rs , g   G k , r; rs , g e

(3.2) V

(3.3)

dk.

Substituting equation (3.3) into equation (3.2), we obtain I (r, r , rs , rg ) 



 



  A r, k s , k g ; rs , rg m r, k s  k g e



 dk dk , s g

i k s  k g r 

(3.4)

Figure 3.2. Upper panel: 2-D SEG/EAGE salt velocity model; Middle panel: Wavenumber domain resolution functions; and Bottom panel: Enlarged resolution functions for selected locations. where





m r, k s  k g  V m  r   e





 dv ,

i k s  k g r 

(3.5)



A r, k s , k g ; rs , rg  2k02G*  k s , r; rs  G  k s , r; rs 



 



G* k g , r; rg G k g , r; rg .

29

(3.6)

In above equations, k s and k g are local transforms with respect to r or r (not rs and rg ), subscripts s and g are for source- and receiver-side wavefields, respectively. In equation (4), introducing transforms k d  k g  k s and k c  k g  k s , integrating A  r, k d , k c , rs , rg  with respect to k c and recovering the frequency variable  , we have I ( , r, r , rs , rg )





 s    A  , r, k d ; rs , rg m  r, k d  eik d r dk d , 

(3.7)

where s   is the source spectrum, A , r, k d , rs , rg  is the wavenumber domain illumination function given by Xie, et al. (2004). For an acquisition system composed of multiple sources and receivers and within a finite frequency band, the image can be expressed as I  r, k d   R  r, k d  m  r, k d  , (3.8) where R  r, k d    s     A  , r, k d , rs , rg d  . (3.9) rs rg

Figure 3.3. Upper panel: Point scatters used to test the resolution function. Lower panel: Point spreading functions in the SEG/EAGE salt model. The center frequency of the Ricker wavelet is 15 Hz.

30

Figure 3.4. Upper panel: Small 40-degree dipping structures used to test the resolution function. Lower panel: The migration images of the dipping structures in the SEG/EAGE salt model. The center frequency of the Ricker wavelet is 15 Hz.

Figure 3.5. The process of generating a synthetic image from the velocity model. Shown in the upper panel is velocity model, and in the lower panel are localized operations from velocity model to synthetic image. The dashed square indicates the region under processing.

31

Figure 3.6. Upper panel: Reflectors in the SEG/EAGE salt velocity model. Lower panel: Synthetic prestack depth image. Equation (8) can be expressed in the space domain as well

I  r, r    R  r, r    m  r, r   ,

(3.10)

where “  ” stands for the convolution with respect to r  and  R (r, r )   R  r, k d  eik d r dk d . (3.11) Equations (3.9) and (3.11) give the point spreading function in wavenumber and space domains. They are related by the Fourier transform. Equation (3.9) also shows the relationship between the illumination and resolution. We see from equations (3.8) and (3.10) that the image I is a distorted version of the model m . These equations show that, given the resolution function R , one can generate a synthetic “migration image” by convolving R with the model m (Lecomte, et al. 2003), or correct the illumination caused image distortion by deconvolving R from the migration images (Sjoeberg and Lecomte, 2003).

3.3. Numerical Examples To demonstrate the resolution analysis using a wave-equation propagator, a group of numerical examples are calculated using the 2-D SEG/EAGE salt model and its acquisition configuration. The wavefield is extrapolated using a wide-angle one-way propagator (Xie and Wu, 1998). Shown in Figure 3.2 are the velocity model and the wavenumber domain resolution function using equation (3.9). The enlarged resolution functions are shown for selected locations in the bottom panel of Figure 3.2. In the sub-salt region, the resolution function is weak and distorted, resulting in poor-quality images for structures oriented in certain directions. Two examples are calculated to show the responses of structures to the background velocity model and acquisition system. The upper panel of Figure 3.3 shows an array of point scatters embedded in the SEG/EAGE salt model. The resulting images in the lower panel of Figure 3.3

32

simulate the point spreading functions at different locations in the model. In Figure 3.4, an array of 40-degree dipping reflectors is used. The migration image reveals the illumination to the specific dipping angle.

Figure 3.7. The process of illumination correction for the prestack depth image. Upper panel: Prestack image of the salt model; Lower panel: localized operations for image correction. The dashed square indicates the region under processing.

Figure 3.8. Upper panel: The original prestack depth image for the salt model. Lower panel: Resolution deconvolved depth image. Figure 3.5 explains the process of generating a synthetic image. The upper panel is the velocity model. The region within the dashed square is chosen to demonstrate the precess. In the lower panel and from left to right, the velocity model is Fourier transformed to the wavenumber domain and filtered by the wavenumber domain resolution function. Then, the filtered model is 33

transformed back to the space domain. Repeating this process for the entire model generates the synthetic image. Shown in Figure 3.6 are the reflectors and synthetic image generated using this technique. Figure 3.7 explains the process of image correction. The upper panel is the prestack depth image. The region within the dashed square is chosen to demonstrate the processing. In the Lower panel and from left to right, the depth image is Fourier transformed to the wavenumber domain and the resolution function is deconvolved from the image. The process is actually conducted in the wavenumber domain by division. Proper water level is used to maintain the stability. The modified spectrum is transformed back to the space domain and forms the modified image. Repeating this process for the entire model generates the resolution deconvolved image. Figure 3.8 shows the original depth image and the resolution deconvolved image. The images of the subsalt structures, particularly some of the steeply dipping structures, are improved.

3.4. Conclusions In this research, we developed a method for resolution and illumination analysis. Using the wave-equation-based propagator and plane wave decomposition, the method can be applied to complex models without smoothing velocity. Different acquisition configurations can be adopted in the calculation. The resolution function can be used to simulate a synthetic migration image, correct the depth image, and conduct illumination analysis for specific acquisition systems and velocity models.

34

PART II. AMPLITUDE FIDELITY OF SEISMIC IMAGING: TRUE-REFLECTION IMAGING

35

4. OPTIMIZATIONS IN TRUE-AMPLITUDE ONE-WAY PROPAGATORS In the traditional one-way propagators, much less attention is paid to the amplitude than to the phase. They cannot provide accurate amplitude even at the level of leading order asymptotic WKBJ or ray-theoretical amplitudes (Zhang, et al., 2003). WKBJ amplitude was then introduced into the original one-way wave propagator (e.g. Clayton and Stolt, 1981). Traditionally WKBJ solution was derived by asymptotic approximation in smoothly varying v(z) media (e.g. Morse and Feshbach, 1953; Aki and Richards, 1980; Clayton and Stolt, 1981; Stolt and Benson, 1986). It has also been obtained by approximately factorizing the full-wave operator into one-way wave operators in heterogeneous media (Zhang, 1993). With an extra amplitude term introduced to the conventional one-way wave equations, the first-order transport equation of the one-way wave equation is the same as that from the full-wave equation in the sense of high-frequency approximation. In this sense, they are called “true-amplitude” one-way wave equations (Zhang, et al., 2003). It was also derived from the conservation of energy flux in smoothly varying v(z) media and extended to general heterogeneous media in local wavenumber(angle) domain by introducing the concepts of “transparent boundary condition” and “transparent propagator” (Wu & Cao, 2005; Cao & Wu, 2005, 2006; Luo et al., 2005). The concept of transparent propagator is similar to the flux-normalized propagator for general heterogeneous media (Wapenaar & Grimbergen, 1996). In the traditional one-way propagators, the square root operator is approximated by a series of rational functions. One approach to achieve better accuracy for the large propagation-angle waves is to use higher-order expansion. However, it will bring extra computation cost, so here we focus on the other approach, to optimize the expansion coefficients in the rational functions, which has been studied extensively in the traditional one-way ‘phase’ propagator (e.g., Lee & Suh, 1985; Ristow & Rühl, 1994; Huang & Fehler, 2000; Xie et al., 2000). For the conventional 45º FD equation Lee & Suh (1985) proposed a least-squares optimization scheme. Ristow & Rühl (1994) proposed a locally optimized scheme for the Fourier Finite-Difference (FFD) propagator. Huang & Fehler (2000) proposed a globally optimized scheme for the Wide-Angle Screen (WAS) propagator or GSP (Generalized Screen Propagator). Considering the advantage of dual-domain extrapolators (e.g., WAS, FFD) compared with pure FD extrapolators, and the fact that the optimized WAS propagator has similar accuracy to the optimized FFD propagator (Ristow & Rühl, 1994; Huang & Fehler, 2000), we use the globally optimized WAS propagator as the one-way extrapolator here. In this paper, three implementations for the amplitude correction in the true-amplitude one-way wave equation are discussed. Then the optimization schemes for every approximated term in TAOP are applied to these implementations to improve the phase and amplitude accuracy for the large propagation-angle waves. Finally numerical results are shown, including the amplitude error caused by phase distortion, comparisons of different implementations for TAOP with optimizations in smoothly varying heterogeneous media, and the application of TAOP in the smoothed 2D SEG/EAGE salt model.

36

4.1. Implementations of the true-amplitude one-way wave equation The true-amplitude one-way wave equation for generally heterogeneous v(x, y, z) media can be written as (Zhang et al., 2003; 2005)   ( z  i  ) P( xT , z ,  )  0 , (4.1)   P( xT , z  0,  )   1  1 ( x  x s )  2i where xT=(x, y), P(xT, z, ) is pressure, and

 (xT , z ) 

2 v 2 (xT , z )



2 , xT2

1 2   2  2  1 v ( x T , z )   2 2  ( x T , z )  1     v  v  , 2v z xT2   x T2      2 2 2     2 2, 2 xT x y where  is the square-root operator for the traditional one-way propagator and  is the differential operator for the WKBJ correction factor.

(4.2) (4.3)

The straightforward way to solve TAOP (4.1) is, at each depth step firstly solve the traditional one-way wave equation without the amplitude correction term    ( z  i ) P( xT , z,  )  0 ; (4.4)   P( xT , z  0,  )   1  1 ( x  x s )  2i then apply the amplitude WKBJ correction,  (4.5) P (xT , z,  )  P( xT , z,  ) . z

To solve the traditional one-way wave propagator (4.4), first-order Padé approximation can be used to expand the square root operator  (Collins, 1989; Xie & Wu, 1998),  a1 2    k 2 (xT , z ) xT2  , (4.6)  (xT , z )  k ( xT , z )  1  2  b1  1   k 2 ( xT , z ) xT2   where, k (xT , z )   v(xT , z ) . (4.7) The original coefficients are a1=0.5 and b1=0.25. The first-order Padé approximation (4.6) is also the equation for the conventional 45º FD migration. By neglecting the fourth order partial derivatives, the square root operator (4.6) can be further approximated to the dual-domain WAS propagator (Xie & Wu, 1998)

37

2

2  v0 ( z )   2 a ( m  1)   2 2     1   2    xT , (4.8) WAS (xT , z )      1     2 2 v0 ( z )  m  v0 ( z ) 2  v0 ( z )  xT 2  v0 ( z )  1  b2 (1  m )   2    xT in which v0=v0(z) is the reference velocity at depth level z and m(x, z)=v(x, z)/v0(z) is the velocity contrast. The coefficients are originally set as a2=a1 and b2=b1 in (4.6). In general case the surface source term correction can be implemented by FD method with the expansion (4.6).

One-step FD WKBJ correction in space domain The direct way to solve the WKBJ correction (4.5) is using space domain FD method, 1

1  v ( xT , z )  v 2 ( xT , z )  2  P( xT , z,  )  P(xT , z,  ) , (4.9) 1   2 2v ( x T , z ) z z xT2   which is similar to the FD method solving the traditional one-way wave eq. (4.4) and the wideangle FD correction in WAS and FFD propagator. Their computation costs should be similar. Two-step WKBJ correction in dual (space and wavenumber) domain The advantages of above one-step FD implementation for WKBJ correction are that no reference velocity is required which makes this implementation method compatible with any traditional one-way extrapolator; and the FD correction method in space domain should be stable. However, considering the disadvantage of pure FD extrapolators compared with the dual-domain extrapolators, an alternative dual domain two-step method can be applied to the WKBJ correction too. The WKBJ correction (4.5) can be implemented by solving following two equations sequentially,  (4.10) P (xT , z,  )  0 P (xT , z,  ) , z  P(xT , z,  )  (   0 ) P(xT , z,  ) , (4.11) z where, 1

1 v0 ( z )  v02 ( z )  2   ln  0 ( z,  ) 1 2 0    2  2v0 ( z ) z   xT  z



1 2

,

(4.12)

1

 2  2 (4.13)  0 ( z )   k0 2 ( z )  2  , x  T   k0 ( z )   v0 ( z ) , (4.14) are defined for the reference v0(z) model. With eq. (4.12), the first step correction (4.10) can be implemented in the wavenumber domain as below, 1

 k 0 ( z,  )  2  P ( k x , k y , z  z,  )   0 z  P(k x , k y , z,  ) ,  k z ( z  z ,  ) 

38

(4.15)

1 2

k ( z,  )   k0 ( z )  k  k  , (4.16) where k z0 ( z,  ) is the global vertical wavenumber in the reference v0(z) media; kx and ky are the horizontal wavenumbers. The second step is similar to the one-step FD correction (4.9), just substituting  with -0. This dual domain two-step WKBJ correction can easily be implemented in the traditional dual domain one-way propagators (e.g., WAS, FFD). 0 z

2

2 x

2 y

Target-oriented flux transparent propagator Above two implementations for the WKBJ correction need do the correction at every depth step to the target area, which will be time-consuming if we are interested in the amplitude behavior of the pressure only in a very deep target, e.g. sub-salt area. By the variable substitution in TAOP (4.1), another form of TAOP can be obtained (Zhang et al., 2005),   ( z  i ) F ( xT , z ,  )  0 , (4.17)  1  1/ 2  F ( xT , z  0,  )     ( x  x s )  2i where F  1/ 2 P . (4.18) The form of TAOP (4.17) is same as the traditional one-way wave eq. (4.4) except with a different source term modification. The amplitude correction term  is not explicitly included in this equation. This corresponds to using a point source of energy-flux instead of pressure source to excite the wavefield. The pressure field P can be recovered with (4.18) by FD method. With this new scheme, the target-oriented pressure field can be obtained with relatively low cost. In the v(z) media, eq. (4.18) will be F=kz1/2P. Since the energy flux incident on a unit area of the horizontal interface is kzP2/, the downward continuation of F will make the energy flux conservative during the propagation. We call this propagator flux transparent propagator (Wu & Cao, 2005), which is the WKBJ solution for smooth media. It should be noted that theoretically there is no any approximation when transforming the original form of TAOP (4.1) to the energy-flux transparent eq. (4.17). However, for the implementation, we need to approximately expand the fourth root operator 1/2 with some rational functions. This introduces further approximation to the WKBJ differential operator (4.3), which will be shown by the following numerical results too. The first-order Padé approximation for the fourth root operator 1/2 can be written as,  a3 2    k 2 ( xT , z ) xT2  , (4.19) 1/ 2  k ( xT , z )  1  b3 2   1   k 2 ( xT , z ) xT2   with a3=0.25, b3=0.375.

4.2. Optimizations for different terms in the true-amplitude one-way propagators

39

Different one-way propagators have different approximations to the true square root operator (4.2): e.g., for the 45º FD the approximated operator is (4.6); for WAS the approximated operator is (4.8). Though eq. (4.8) is derived from (4.6), a different optimization should be applied for (4.8) since further approximation was applied in getting (4.8). The error may be much bigger than expected for large propagation-angle waves if the optimized coefficients for (4.6) are applied to (4.8) (Figure 4.1). Compared with traditional one-way propagators in which only the phase term needs optimization, it is more complicated for TAOP, in which every approximated term needs optimization. Furthermore, we may have differently approximated operators for the traditional one-way (phase) propagation term, source modification term and the amplitude WKBJ correction term. For example, for the flux transparent propagator implemented with WAS propagator, the WAS propagator is approximated with (4.8); while the source modification term and the pressure recovering term are approximated with (4.19). Different optimization scheme can be applied to each approximated term in TAOP. In the following numerical results, for the WAS approximation WAS, eq. (4.8), the expansion coefficients obtained from the globally optimized scheme by Huang & Fehler (2000) are used; for the first-order Padé approximation to 1/2, eq. (4.19), we use a similar scheme to Huang & Fehler’s (2000) to obtain the optimized coefficients; while for the first-order Padé approximation to , eq. (4.6), we use the expansion coefficients derived from the least-squares optimization scheme by Lee & Suh (1985). After the optimization, the relative errors for above 3 approximated operators are much smaller for the larger propagation-angle waves (Figure 4.2).

relative error (%)

. Figure 4.1: Comparison of the relative errors for differently approximated kz with optimization: The red line is for the first-order Padé approximation to kz (4.6) with its optimized coefficients by Lee & Suh (1985); the green line is for WAS approximation to kz (4.8) with its optimized coefficients by Huang & Fehler (2000) (for maximum perturbation 2); the blue line is for approximation (4.8) but with optimized coefficients for (4.6) by Lee & Suh (1985). 10 9 8 7 6 5 4 3 2 1 0 -1 0

Pade to kz (Original) Pade to kz (Lee & Suh) WAS (Original) WAS (Huang & Fehler) Pade to sqrt(kz) (Original) Pade to sqrt(kz) (Optimized)

10

20

30

40

50

o

60

70

propagation angle ( ) Figure 4.2: Comparison of the relative errors before (dashed lines) and after (solid lines) optimization for three approximated operators: first-order Padé approximation to the square root

40

operator (4.6) (red lines), WAS approximation to the square root operator (4.8) (green lines), and first-order Padé approximation to the fourth root operator (4.19) (blue lines).

4.3. Numerical results For simplicity, the numerical examples shown here are in 2D models, but it is very easy to extend them to 3D. The model dimension is 10.24 km x 5.12 km. The source is located at xs=5.12 km on the surface. The source time function is Ricker wavelet with peak frequency 15 Hz. Firstly, a homogeneous model is used to study the influence of phase distortion on the amplitude of the wavefield in the traditional one-way propagator. Then above discussed three optimized implementations, one-step FD method, two-step dual domain method and flux transparent method, are compared in a smoothly varying heterogeneous model. Finally the WKBJ-corrected amplitude is investigated in the smoothed 2D SEG/EAGE salt model. Phase distortion caused amplitude error In the traditional one-way propagator, much less attention is paid to the amplitude than to the phase. Before the amplitude WKBJ correction is applied to the traditional one-way propagator, we first investigate the amplitude behavior related to the phase errors in the traditional one-way propagator. A homogeneous 4.5 km/s velocity model is used and the velocity contrast m is chosen as 1.2. Figure 4.3 shows the phase and amplitude behavior for phase-screen and WAS propagator compared with the true values. For phase-screen propagator (Figure 4.3a, b), the phase for 0º propagation-angle wave is correct, but its error increases with the propagation angle bigger than 0º. The amplitude is incorrect even for the 0º propagation angle. It is caused by the phase error of waves with propagation-angle bigger than 0º, which distorts the energy flux of the 0º propagation-angle wave. With the wide-angle phase more accurate WAS propagator, the phase of 60º propagation-angle wave is correct and there is a small phase error for 75º (Figure 4.3c). For the amplitude, it is still correct for 45º propagation angle, but a small error occurs for 60º and the error is much bigger for 75º (Figure 4.3d). Above numerical results demonstrate that for the traditional one-way propagator the maximum propagation angle at which the amplitude is correct is smaller than that at which the phase is correct. Therefore, to achieve better amplitude accuracy for large propagation-angle waves, optimization should also be applied to the traditional one-way phase propagator. Figure 4.4 shows the results for a larger velocity contrast m=1.62. With this larger contrast, the phase errors from both phase-screen and WAS propagator increase for the wave with large propagation angle, which makes the amplitude error caused by the phase distortion increase too. Comparison of different implementations for optimized true-amplitude propagators In this part, three implementations with optimization are compared in a velocity model varying linearly in the vertical direction but as a sine function in the horizontal direction: v(x, z)=3.0+0.18*sin(2π(x/xs-1/4))+0.36z (km/s). The extracted amplitudes along the radial direction 0º to 75º for different schemes are shown in Figure 4.5. The amplitude before the WKBJ

41

correction has obvious error compared with the true amplitude which is obtained from the fullwave FD method (Figure 4.5a). With the one-step FD correction in space domain the amplitude is corrected very well up to 60º (Figure 4.5b). With the two-step method the amplitude is corrected very well up to 75º (Figure 4.5c), which is better than the one-step method. Figure 4.5d and 5e are for the flux transparent implementation. The corrected pressure is better after the optimization to the fourth root operator (2.19). However, even with optimization, the corrected amplitude for 75º is still not as good as the other two methods based on the original differential operator (2.3). This is consistent with the theoretical analysis. Application in smoothed 2D SEG/EAGE salt model Finally, the WKBJ-corrected amplitude is investigated in a more realistic model, 2D SEG/EAGE salt model. Here the salt model is smoothed to make it transparent to the waves (Figure 4.6a). A group of receivers are aligned horizontally at the bottom in the sub-salt area (green line in Figure 4.6a). The amplitudes from one-step method and transparent method both agree very well with the true amplitude (Figure 4.6b). However, there is almost no extra cost for the transparent method to obtain above WKBJ-corrected amplitude since only the pressure at the receivers’ depth level needs to be recovered from the downward-continued flux field.

42

Figure 4.3: Phase and amplitude behavior for the phase-screen and WAS propagator with a constant velocity contrast m=1.2. (a), (b) are for phase-screen and (c), (d) are for WAS. (a, c): impulse responses at 3 time steps; the true wavefronts are shown with red lines; and the green lines show the receivers along radial directions 0º to 75º with a 15º interval. (b, d): picked maximum amplitude versus distance from the source for radial directions shown in (a, c); the solid lines are the results from phaseshift propagator with true velocity and the dashed lines are results from phase-screen or WAS propagator.

Figure 4.4: Same as Figure 4.3 but with a larger velocity contrast m=1.62.

43

D No WKBJ, GOWAS, Pade with Lee&Suh’s a&b

E kwkbj=11, GOWAS, Pade with Lee&Suh’s a&b

( v=3+0.18sin(2π(x/5.12-1/4)+0.36z (km/s )) 75 60 45 30 15 0

Absolute amplitude

0.03 0.025

0.1

( v=3+0.18sin(2π(x/5.12-1/4)+0.36z (km/s ))

0.02

0.15

0.015 1

0.05

2

3

4

Solid lines: Full wave FD Dashed lines: One way without WKBJ

0 0

1 2 3 Distance from the source ( km )

0.1

1

0.025

0.1

0.02

0.15

0.015 1

0.05

0 0

2

3

4

Solid lines: Full wave FD Dashed lines: One way with WKBJ

0.15

1 2 3 Distance from the source ( km )

0.035

75 60 45 30 15 0

0.025

0.1

0.02 0.015 1

2

0.05

Absolute amplitude

3

4

Solid lines: Full wave FD Dashed lines: One way with WKBJ

0 0

4

4

kf2p=1, GOWAS, Pade with original a&b ( v=3+0.18sin(2π(x/5.12-1/4)+0.36z (km/s ))

0.035

4

75 60 45 30 15 0

0.03 0.025

0.1

4

1 2 3 Distance from the source ( km ) kf2p=1, GOWAS, GO for sqrt(kz) ( v=3+0.18sin(2π(x/5.12-1/4)+0.36z (km/s ))

1 2 3 Distance from the source ( km )

H

3

Solid lines: Full wave FD Dashed lines: One way with WKBJ

0.03

Absolute amplitude

Absolute amplitude

0.03

2

0.05

G

75 60 45 30 15 0

0.02 0.015

( v=3+0.18sin(2π(x/5.12-1/4)+0.36z (km/s )) 0.035

75 60 45 30 15 0

0.025

0 0

4

F kwkbj=21, GOWAS, Pade with Lee&Suh’s a&b 0.15

0.035 0.03

Absolute amplitude

0.15

0.035

0.02 0.015 1

0.05

0 0

2

3

4

Solid lines: Full wave FD Dashed lines: One way with WKBJ

1 2 3 Distance from the source ( km )

4

Figure4.5: Comparison of the WKBJ correction on the amplitudes for the three implementations in a smoothly varying medium v(x, z) =3.0+0.18*sin(2π(x/xs-1/4))+0.36z (km/s). Different colored lines are for different radial directions =0º, 15º, 30º, 45º, 60º, 75º. The solid lines are the results from full-wave FD method, and the dashed lines are results from one-way WAS method. The inset figures are zoomed in versions of the main figures. (a) No WKBJ; (b) WKBJ with onestep FD correction; (c) WKBJ with two-step correction; (d) flux transparent implementation without optimization to (2.19); (e) flux transparent implementation with optimization to (2.19).

44

(a)

Absolute amplitude

Effect of WKBJ on amplitude, GOWAS (smoothed SEG/EAGE salt model) 0.012 0.01 0.008 0.006 0.004 0.002 0 0

Full−wave FD One−way without WKBJ One−step WKBJ Transparent propagator

1 2 3 Distance from the center ( km )

(b) Figure 4.6: Amplitude behavior with/out the WKBJ correction for the smoothed 2D SEG/EAGE salt model: (a) velocity model with the source at the surface (star) and some receivers aligned at the bottom in sub-salt area (green line); (b) picked maximum amplitude with/out the WKBJ correction along the receivers shown in (a).

4.4. Conclusions Three implementations (one-step FD method, two-step dual domain method and target-oriented flux transparent method) for the amplitude correction in TAOP are investigated. In the traditional one-way propagator, the phase distortion makes the maximum propagation angle at which the amplitude is correct smaller than that at which the phase is correct. To achieve better wavefield accuracy for the large propagation-angle waves, optimization schemes are applied to all the approximated terms in TAOP including the traditional one-way (phase) propagation term, surface source modification term and also the amplitude WKBJ correction term for the flux transparent method. What makes it more complicated is that each term may be approximated differently so that different optimizations should be applied correspondingly. Theoretical analysis and numerical results in a heterogeneous model show that the three implementations are different in accuracy, efficiency and compatibility with traditional one-way extrapolators: twostep method is more accurate than one-step method but one-step method is compatible with any traditional one-way extrapolator; the flux transparent method is less accurate than one-step method due to the approximated implementation to 1/2; after optimizing the expansion for 1/2 45

its accuracy is much more close to the one-step method; the advantage of the flux transparent method is that target-oriented analysis can be applied with it, which may save a lot of computation. In the smoothed 2D SEG/EAGE salt model, the amplitude distribution beneath the salt dome given by TAOP agrees well with the true amplitude from full-wave FD method, and there is almost no extra cost for the flux transparent method to obtain this WKBJ-corrected amplitude for that target.

46

5. SUPERWIDE-ANGLE ONE-WAY MIGRATION: APPLICATION TO 2004 BP 2D BENCHMARK MODEL 5.1. Background The one-way wave equation simulates the exact one-way wavefield only in very limited propagation directions. Quite a few methods have been proposed to improve the wide-angle accuracy of one-way wave propagators (Collins and Westwood, 1991; Ristow and Ruhl, 1994; Wu, 1994, 1996; Xie and Wu, 1998; Grimbergen et al., 1998; Jin et al., 1998, 2002; Huang et al., 1999a, b; Xie et al., 2000; De Hoop et al., 2000; Le Rousseau et al., 2001; Han and Wu, 2005; Thomson, 2005). In spite of these progresses, the accuracy of wide-angle waves for strong contrast media is still a serious problem and put a practical limitation on applying these methods to steep reflector imaging. The superwide-angle one-way method (Wu and Jia, 2006) has been developed to improve the regular one-way propagators by a wavefront reconstruction scheme which combines and interpolates the weighted orthogonally propagated one-way wavefields to rebuild the accurate wavefront. The scheme employs the wavefield gradients to determine the weighting function for the two wavefields. In this paper, we will show the capability of this method on imaging the steep and overhanging salt flanks through the example of 2004 BP 2D benchmark model dataset.

5.2. Theory Assume two one-way wave equations in z- and x-directions as u D  x , z   PD u D  x, z  , z u H  x , z   PH u H  x , z  , x where uD  x, z  and uH  x, z  are the one-way wavefields propagated downward in z-

(5.1)

direction and horizontally in x-direction, respectively; PD and PH are the corresponding one-way propagators. The solution for one step forward propagation with respect to a reference point x0   x0 , z0  can be written as uD  x, z0  z   PD uD  x, z0   , uH  x0  x, z   PH uH  x0 , z   ,

(5.2)

where x and z are step sizes, incorporated into the propagators on the right-hand side of eq. (5.2). Now considering the one-way propagations in z-direction having n steps forward and in xdirection having m steps forward, we have uD  x, z0  nz   PDn uD  x, z0   , uH  x0  mx, z   PHm uH  x0 , z   .

47

(5.3)

We reconstruct the accurate wavefield by taking a weighted average of the two one-way wavefields. Wavefront reconstruction is a stepwise marching process, which means the wavefront reconstructed at the current step will be used for further one-way propagation. Note that the one-way propagation is exact in the forward direction and very accurate for small-angle waves. Therefore, at any location we put heavy weight on the one-way solution which propagates small-angle waves at that point and light weight on the other one-way solution which propagates large-angle waves. A weighting scheme is constructed according to the local propagating angles by (5.4) u  x, z   wD  D  uD  x, z   wH  H  uH  x, z  , where u  x, z  is the reconstructed wavefield,  D is the propagating angle with respect to the zaxis, and  H , to the x-axis; wD  D  and wH  H  are the weights for uD  x, z  and uH  x, z  , respectively. Eqs (5.3) and (5.4) constitute the complete scheme of wavefront reconstruction. The numbers n and m in eq. (5.3) determine how often wavefront reconstruction is conducted. If n=m=1, the interpolation shown in eq. (5.4) will be implemented immediately after each step of the two one-way propagations. In other words, it is the case of stepwise continuous reconstruction. On the other hand, we can give n=N and m=M in eq. (5.3), where N and M are the total sampling numbers of the propagation area in z- and x-directions, respectively. This corresponds to a simple scheme of weighted average without updating the wavefront interactively. The interpolation of the two one-way wavefields is separated completely from the wave propagations. As a result, the propagations remain the same as the regular one-way wavefield propagations. Given n=N and m=1, we have an intermediate case of the reconstruction scheme. In this case, the downward wavefield is propagated conventionally; then at each step of the horizontal propagation, the wavefield is updated by weighted summation of the two corresponding wavefields, in the way of eq. (5.4). We employ this intermediate scheme to reconstruct the wavefield in the following numerical experiments. The weight function is a key point to the wavefront reconstruction method. The wavefield gradients are used to calculate the weights for the two one-way wavefields. We can use either horizontal screens (i.e. downward one-way wave propagation) or vertical screens (i.e. horizontal one-way wave propagation) to obtain the wavefield gradients. Since neither of these two oneway propagators is accurate at large propagation angles, the final wavefield gradient at a given point should be determined by weighted summation of the gradients obtained from the two orthogonal one-way wave propagations.

5.3. Experiments on 2004 BP 2D Benchmark Dataset The weighting summation can be applied to most of the one-way propagators. In this work, we test the scheme using the wide-angle Padé GSP algorithm (Xie and Wu, 1998) and the LCB beamlet propagator (Wu et al., 2000; Wang and Wu, 2002; Luo et al., 2004), respectively. In this section the superwide-angle one-way migration will be demonstrated on the 2004 BP 2D dataset (Billette and Brandsberg-Dahl, 2005). Figure 5.1 shows the velocity model under consideration. Our interest is imaging the boundaries of the salt body which has a constant

48

velocity of 4790m/s. The salt has extremely steep and even overhanging dips at its left flanks, which present a great challenge to most one-way migration methods.

Figure 5.1: velocity model of 2004 BP 2D benchmark dataset.

Figure 5.2: modeling results of the superwide-angle GSP method. The source is located at (28 km, 0). Times for these snapshots are 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5 and 6.0 s. The contour of the salt body is sketched. Figure 5.2 shows the snapshots of nine times using the superwide-angle GSP. The shot is located at the surface and the dominant frequency is 27 Hz. Note the phases concerning with the reflected turning waves near the left flank of the salt dome. The synthetic data for imaging were generated by finite difference method with receiver interval of 12.5m, shot interval of 50m and receiver array length of 15km. 14 seconds of data were recorded and the dominant frequency was 27Hz. The signals of turning waves reflected on the salt flank could be traced only in a limited part of the whole single shot gathers. The shot locations range from x=25km to x=30km for these relevant gathers. This can be seen clearly in Figure 5.3 which shows the rays tracing from the salt boundary. The geometry under consideration can record the primary turning reflections on the vertical part of the salt flank. However, for the overhanging flank, no primary reflection can be acquired. Instead, the duplex turning reflections are recorded by some shot gathers with the shot locations near x=28km. Based on these data, Figure 5.5 shows the stacked images using both the regular GSP method and the superwide-angle GSP method. Figure 5.4 shows the single-shot images of the two

49

methods. Figures 5.6 and 5.7 show the corresponding results for the LCB propagator. Due to the angle limitation, the regular downward propagator can hardly simulate turning waves and therefore contributes nothing to the image of the overhanging flank. Note the ‘gap’ between the depths of z=5km and z=8km in the image of the salt body when using regular one-way method for migration. However, in the imaging results of the superwide-angle methods, the phases of the overhanging flank are sharp and clear and their locations are accurate. Due to the local perturbation assumption, the LCB beamlet propagator can not handle the reflections very well. This leads to relatively weak images of the salt flank using the superwideangle LCB method, compared with the results of the superwide-angle GSP method.

Figure 5.3: rays normally shooting from the salt boundary. The contour of the salt body is sketched partly. (Provided by Yingcai Zheng). If we implement an extra downward one-way pass in addition to the regular superwide-angle scheme, the image will be further improved. Figure 5.8 shows the image obtained in this way. We use the data with the shot locations from x=25km to x=40.8km. The image of the salt boundary is more continuous. The energy has been rebuilt on almost all boundaries of the salt body.

5.4. Discussions and Conclusions Although the wavefront reconstruction method yields quite good images of the steep and overhanging features, it can still be improved when full reconstruction scheme is implemented. Full reconstruction is necessary to reduce the artifacts and simulate the turning waves accurately. The key point of full reconstruction is the implementation of weighting the two one-way wavefields. Wavefield gradients can still be used to determine the weights in this case. When the superwide-angle one-way method is applied to local propagators, such as the LCB beamlet propagator, the weights can also be calculated in the local angle domain, which provides an easy approach to handle the full reconstruction. This work demonstrates the accuracy of the superwide-angle one-way propagator by using wavefield gradients for the calculation of weight functions. The method is also applied to seismic imaging for steep salt flanks using turning waves. Numerical examples show that the superwideangle one-way method can model large angle waves well beyond 90º. The image quality, 50

especially for overhanging salt flanks, has been markedly improved compared with the standard one-way propagation methods.

(a)

(b) Figure 5.4: comparison of the single-shot images obtained by the regular GSP method and the superwide-angle GSP method with the shot location at x=28km. (a) regular downward GSP method; (b) superwide-angle GSP method.

51

(a)

(b) Figure 5.5: comparison of the stacked images obtained by the regular GSP method and the superwide-angle GSP method with the shots from x=25km to x=30km. (a) regular downward GSP method; (b) superwide-angle GSP method.

52

(a)

(b) Figure 5.6: Comparison of the single-shot images obtained by the regular LCB method and the superwide-angle LCB method with the shot location at x=28km. (a) regular downward LCB method; (b) superwide-angle LCB method.

53

(a)

(b) Figure 5.7: comparison of the stacked images obtained by the regular LCB method and the superwide-angle LCB method with the shots from x=25km to x=30km. (a) regular downward LCB method; (b) superwide-angle LCB method.

54

Figure 5.8: image obtained by superwide-angle GSP plus downward GSP method.

55

6. TARGET ORIENTED FULL-WAVE EQUATION BASED ILLUMINATION ANALYSIS 6.1 Background The seismic illumination analysis quantifies the capability that an acquisition system can image the subsurface structures. The illumination measurements can be used to optimize the acquisition survey design, evaluate the image quality, and make corrections to seismic image, resulting in more accurate subsurface physical parameter retrieval. Although the early attempts in illumination estimates adopted oversimplified model, e.g., homogeneous velocity model, horizontal target and symmetric ray paths, recently it has developed quickly. The illumination calculation can be considered a special type of modeling. First, seismic waves from sources and receivers need to be propagated through overburden structures to reach to the target. Second, the propagation directions for both waves need be determined at the target location since targets with different dipping angles can only be illuminated by waves at certain directions. Finally, the illumination measurement can be calculated from incident-scattering waves from different directions. The ray-based methods can handle smoothly varying heterogeneous velocity models and has no angle limitations. Particularly, because the ray method naturally carries the wave propagation directions, it has been widely used for illumination calculations (Schneider and Winbow, 1999; Muerdter and Ratcliff, 2001; Gelius, et al., 2002). However, the high-frequency asymptotic approximation causes intrinsic difficulties for this method being used in complex overburden structures. Recently, the one-way propagator has been widely used in seismic imaging. At the same time, methods specifically for extracting angle information from the wavefield are developed (Xie et al., 2002; Wu and Chen, 2002, 2003; Sava and Fomel, 2003). Combining the two techniques together, a one-way wave equation based method is developed for illumination analysis (Xie et al. 2003, 2005, 2006; Wu and Chen 2002, 2006; Wu et al. 2003; Jin et al. 2003, 2006; Mao and Wu, 2007; Cao and Wu, 2008a). This is an efficient approach and is consistent with most oneway wave equation based migration methods. However this approach does not properly handle the wide-angle signals such as turning waves. In the last few years, reverse-time migration regains the attention in seismic community partially because it does not have the angle limitation and can properly image the steeply dipping structures with turning waves. It is also a challenge to conduct the illumination analysis that is consistent with the reverse-time migration. To overcome the angle limitation of the one-way propagator, Xie and Yang (2008) proposed a method which uses the full-wave finite-difference method as the propagator and uses a time-domain local slowness analysis method to determine the angle information and calculate the illumination. Cao and Wu (2008b) proposed another fullwave method which uses frequency domain split-step local plan-wave decomposition for angle analysis. In this paper, we further develop the method by Xie and Yang (2008) and present the target oriented illumination analysis using the full-wave propagator. The new method does not 56

have angle limitations usually encountered by a one-way approach. In addition, it accurately solves the energy transmission in a complex velocity model. This method is particularly useful in providing illumination analysis for reverse-time migration and VSP data imaging.

Figure 6.1: Cartoon illustrating the geometry of the illumination analysis.

6.2 Formulations for Different Illumination Measurements The illumination analysis for a target can be started from investigating an incident source beam towards the target at angle  s and a scattered beam leaving the target towards the receiver at angle  g . This pair of incidence-scattering beams forms a basic illumination measurement (see Figure 6.1). Note, due to multi-pathing, scattering and reflection, the waves from a source can come to a target from multiple directions, and so is the case for waves leaving the target back to the receiver. We collect all possible incidence and scattering beams from one pair of sourcereceiver and they give complete description of illumination for the target. In frequency domain, such collection of beams can be expressed with a local illumination matrix (Xie, et al. 2006)





A , r,  s ,  g ; rs , rg 

2 2 v2



I ,  s , r; rs  I ,  g , r; rg



(6.1)

where  is the frequency, v is local velocity, I is the energy flux in  direction, r is the target location, rs and rg are source and receiver locations respectively. Equivalently, we can calculate a similar local illumination matrix using broad band time domain wavefield. Since the timedomain wavefield does not directly provide the propagation direction, we first decompose the Green’s function into angle domain beams using a method similar to that used by Xie et al. (2005) (6.2) G  s , r; rs , t   C  W  r  r  G r ; rs , t  p s   r   r   dr 

where G  r; rs , t  is the Green’s function generated by the FD method, G  s , r; rs , t  is its angledomain decomposition, p is the slowness vector, W is a space window centered at r , C is a normalization factor. A similar equation can be found for scattering beam G  g , r; rg , t  . Considering the frequency and time domain equivalents

57

2

 2 I ,  s , r; rs     t G  s , r; rs , t   dt









I ,  g , r; rg   G 2  g , r; rg , t dt

,

(6.3)

the local illumination matrix can be expressed as



 

A r ,  s ,  g ; rs , rg  2 v 2

  

 t

2





G  s , r; rs , t   dt   G 2  g , r; rg , t dt .

(6.4)

For an acquisition system composed of multiple sources and receivers, its local illumination matrix can be obtained by summing up the contributions from individual source-receiver pairs A  r ,  s ,  g     A  r ,  s ,  g ; rs , rg  . (6.5) rs rg

where Ar,  s ,  g  is the local illumination matrix for the entire acquisition system, and the summation is based on the acquisition geometry, e.g., source and receiver locations, navigation directions, cable lengths. In above equations, the angles  s and  g are related to the acquisition

coordinate. With coordinate transforms  d   g   s  2 and  r   g   s  2 , we can obtain the

illumination matrix A  r ,  d ,  r  defined in target coordinate system, where  d and  r are target dipping angle and reflection angle (see Figure 6.1), respectively. Comparing to the acquisition coordinate, the target coordinate is suitable for investigating the target oriented illumination.

6.3. Numerical Examples for Illumination Calculations We use numerical examples to show the calculations of different illuminations. Local illumination matrix: Figure 6.2 is a sketch illustrating the basic concept of the illumination analysis. Shown in Figure 6.2a is a -45 degree acquisition dip response (ADR) map, calculated in a constant velocity model, overlapped with the source and receiver waves. The slowness analyses shown in 6.2b and 6.2c give the wave propagation directions. At the target location, we construct a local illumination matrix (LIM) (shown in 6.2d) in which the horizontal and vertical coordinates are incidence and scattering angles while the two diagonals are dipping angle and reflection angle, respectively.

Shown in Figure 6.3 are local illumination matrices located at selected locations in a constant velocity model. The sources cover the surface with left sided receiver geometry from 0 to 5.4 km. For a constant velocity model and surface acquisition, the incidence and scattering angles,  s and  g , are limited within 90o . Only the inner square in these illumination matrixes can be illuminated and the angle coverage decreases with increasing depth. The maximum dipping angle which can be illuminated is 90 degrees. Shown in Figure 6.4 are similar illumination matrices in a model composed of a depth dependent background and a salt dome with steep and over-hang flanks. The acquisition system is the same as that used in Figure 6.3. Most of these local illumination matrices cover dipping angles beyond 90o . The local illumination matrices become complex closing to the salt body, partially affected by the reflection waves generated on the salt body. Note that the area can be handled by a full-wave illumination method (the entire outer square) is much larger than that can be handled by a one-way propagator based illumination method (the inner square). This makes the full-wave method can investigate structures with dipping angles beyond 90o (e.g., overhung structures). 58

Figure 6.5 is similar to Figure 6.4 except the illumination matrixes are calculated using most energetic phases picked from time domain. Comparing to Figure 6.4, the apparent difference is there lack of many large angle illuminations from secondary phases.

Figure 6.2: Sketch showing the calculation of illumination measurements, with (a) -45 degree ADR map in a constant velocity model and from a source-receiver pair, (b) source-side and (c) receiver-side slowness analyses, (d) local illumination matrix and (e) illumination as a function of reflection angle for a dipping reflector.

Figure 6.3: Local illumination matrices at selected locations in a constant velocity model generated by multiple shots and receivers. 59

Figure 6.4: Full wave local illumination matrices at selected locations in salt dome model.

Figure 6.5: Most energetic phases local illumination matrices at selected locations in salt dome model.

60

Figure 6.6: Acquisition dip response and local illumination matrix for different source-receiver configuration. The acquisition dip response: The acquisition dip response (ADR) is the total illumination for a target with a given dipping angle. The ADR can be obtained from the illumination matrix by substituting the dipping angle d with the target dipping angle and integrate the matrix over reflection angle  r (Xie and Yang, 2008) D(r,  d )   A  r,  d ,  r  d r . (6.6)

To understand the ADR calculation, a 45o target is added in Figure 6.1a and the 45o ADR is generated by a pair of source-receiver in a constant velocity model. Illustrated in Figure 6.6 are ADR maps and their relations with the illumination matrixes in different velocity models and for source-receiver pairs with different acquisition geometries. The locations of sources and receivers and the sketch of ray paths are illustrated to demonstrate how these factors affecting the illumination. Figure 6.6a shows the -45 degree ADR from a surface acquisition system in a constant velocity model. This configuration cannot properly image a steep structure. Figure 6.6b shows the -120 degree ADR for a surface acquisition in a v(z) velocity model. Because turning waves are generated in this case, a surface acquisition system can image steep or even overhung structures. Figure 6.6c shows the -90 degree ADR for VSP acquisition geometry in a constant velocity model. A VSP acquisition can image a vertical target in a constant velocity model. Shown in the right column in Figure 6.6 are illumination matrixes corresponding to 3 target locations shown in left column. The illumination matrix is calculated using equation (6.4). The energy is located at dipping angle -45 degrees in the right column in Figure 6.6a. In the right column in Figure 6.6b, the dipping angle is located at about -120 degrees which exceeds the 90 degree limitation of the one-way method. Shown in the right column in Figure 6.6c is the illumination matrix for VSP geometry. The energy is located at dipping angle -

61

90 degrees but the reflection angle is slightly away from 0 degree due to the larger reflection angle involved. Shown in Figures 6.7a-d are ADRs for -135, -90, -45, and 0 degrees in a v  z  model and generated by a pair of source-receiver. The model is 15 km in horizontal by 3.6 km in depth and the velocities at the top and bottom are 1.5 and 3.5 km/s, respectively. A 20-Hz Ricker wavelet is used in the FD calculation. As can be seen from the figure, turning waves in a v  z  model illuminate structures with their dipping angle exceeding 90o , which is the limitation of a oneway wave equation based method. To help clearly seeing the illumination at large depths, the automatic gain control (AGC) is used in calculating the response of a single source-receiver pair.

Figure 6.7, ADR maps generated by a source-receiver pair in a v(z) velocity model. Different panels are for different dipping angles and AGC is used for both source and receiver wavefields.

62

Figure 6.8: ADR maps from a source-receiver pair in a two-layer model: (a) -45º ADR for down-down waves; (b) -120º ADR for up-up waves; and (c) -120º ADR for down-up waves.

Figure 6.9: Comparison of ADR maps from different kind of waves by one source-receiver pair in V(z) model and V(z) with a constant basement model: (a) -120º ADR for turning waves in V(z)

63

model; Panels (b)-(d) are in V(z) plus a constant basement from the interface at 3.12 km. (b) 120º ADR for down-down waves; and (c) -120º ADR for down-up waves. (d) -120º ADR for full waves. The FD calculated Green’s functions compose of direct waves, reflections, turning waves and multiples, etc. The full-wave illumination analysis includes contributions from all these waves, while a one-way propagator usually cannot handle up-going waves. Shown in Figure 6.8 are ADRs for a two-layer velocity model. The velocities in the top and bottom layers are 2.5 and 4.0 km/s, and the interface is located at depth 3.12 km. We use this model to analyze contributions from different waves which can be separated by propagation directions and time windows. Shown in Figure 6.8a is the 45o ADR which is generated by direct “down-down” waves. Figure 6.8b is the 120o ADR generated by “up-up” reflections. Figure 6.8c is the 120o ADR obtained from the direct source wave and reflected receiver wave, i.e. the “down-up” waves. The later two ADRs cannot be properly handled by a one-way propagator based illumination analysis. Figure 6.9 compares -120º ADR maps generated by different waves from a single sourcereceiver pair. Figure 6.9a is the -120º ADR by turning waves similar with Figure 6.7 in a V(z) model. Figures 9b-d are -120º ADRs in a model which has a V(z) velocity in the upper layer and a constant velocity below depth 3.12 km. We use a time window to pick the most energetic phases. This method eliminates most of the reflection waves and head waves, except the time interval that those waves interfere together. Those interference wavefronts cause the higher illumination intensity in the bottom part of Figure 6.9b-d. Figure 6.9b is generated by most energetic phases both in source and receiver side waves. In figure 6.9c, the source wavefield contains all waves and the receiver side includes only the most energetic phases. Figure 6.9d is generated by full waves. Notes the illumination area becomes larger from 6.9a-d. Illumination as a function of target reflection angle: For a locally planar target with a normal vector pointing to direction  n , the target-oriented illumination as a function of reflection angle can be obtained from the illumination matrix (Xie, et al., 2006) . (6.7) D  r,  r ,  n   A  r ,  d ,  r    d

n

Figure 6.2d illustrates the calculation of D  r, r ,  n  . The dashed line passes through the energy peak is  d  -45º and the coordinate along the line is reflection angle  r . The illumination as a function of reflection angle can be directly measured in the illumination matrix along this line. As an example, the illumination distribution for a -45º target is shown in Figure 6.2e as a fanshaped distribution. Note the maximum illumination is slightly biased from the normal direction, resulting from the finite offset between the source and receiver (see 6.2a and 6.2d). Figure 6.10 compares the target oriented illumination with the reverse-time image for a salt dome model. Figure 6.10a shows the reverse-time image (Xie and Wu, 2006) where only one shot is used and the receivers cover trace numbers 0 to 320. Shown in 6.10b is the illumination as a function of reflection angles along the bottom reflector and salt flank. The illumination estimate is consistent with the migration image. The horizontal interface and vertical salt flank are properly illuminated and imaged. However, the overhung part of the salt flank (indicated by the arrow) is missing from the image. The target oriented illumination properly indicates that there is an illumination hole.

64

Figure 6.10: Comparison between reverse-time image and target illumination for the salt dome model with (a) reverse-time image from one shot and (b) illumination as a function of reflection angle along the target.

Figure 6.11: ADRs from full waves for different dipping angles in the salt dome model.

65

Figure 6.12: ADRs from most energetic phases for different dipping angles in the salt dome model. The volumetric ADR: Figure 6.11 shows the volumetric ADR for the salt dome model shown in Figure 6.4. We use 46 surface shots to generate the illumination. Figure 6.11a is the -135º ADR map. The overhung structure of left salt flank can be illuminated by the adopted acquisition geometry. The weak illumination from the reflection waves can be found at the left side of the dome structure, and there is a shadow area on the right side of the dome. Similarly, Figure 6.11b indicates that the -90 degree left salt flank can be well illuminated. The shadow area is much smaller than that in Figure 6.11a. Figures 6.11c-d show that structures with small dipping angles can be evenly illuminated. Figure 6.12 is similar to Figure 6.11 except the former is generated by the most energetic arrivals. Comparing with Figure 6.11, there are apparent shadow zones near the salt dome.

Shown in Figure 6.13 is part of the BP velocity model (Billette and Brandsberg-Dahl, 2005) which is composed of a complex background and salt bodies with overhangs. The sources cover the surface of the model. The left sided receivers have offset between 0 to -15 km. Figures 6.14ad are ADRs for 0, -45, -90 and -135 degrees, respectively. The 0 degree ADR in Figure 6.14a is generally strong throughout the model and is consistent with the migration result where the horizontal structures are well imaged. The -90 degree ADR in Figure 6.14c is responsible for illuminating the vertical wall of the salt dome. The -135 degree ADR shown in Figure 6.14d contributes to the illumination of the overhung part of the salt flank. Imaging this overhung structure is a major challenger for one-way wave equation based method, while the reverse time migration can image this structure satisfactorily. The illumination method based on the full-wave 66

propagator is consistent with the reverse time migration and can provide corresponding illumination analysis.

Figure 13: Salt dome structure in 2D BP model.

Figure 6.14: ADRs for different dipping angles in the BP salt dome structure.

6.4. Conclusions We propose a full-wave equation based illumination analysis method which has the target oriented capability. The full-wave finite difference propagator is used to extrapolate the source and receiver waves to the subsurface target. A time domain local slowness analysis is used to decompose the wave fields into local angle beams and form the local illumination matrix. ADRs beyond 90 degrees can be extracted from the illumination matrix. The current approach has no angle limitations and is particularly useful for providing illumination information for reversetime migration.

67

7. LOCAL-ANGLE DOMAIN ILLUMINATION IN FREQUENCY DOMAIN FOR FULL-WAVE PROPAGATORS The illumination of seismic waves to a subsurface target is affected by many factors, e.g., the limited acquisition geometry, the complex overburden structure, and the reflector dip angle. Illumination analysis can be used for survey design (e.g., Li and Dong, 2006). An uneven illumination also produces a distorted image. Illumination analysis is a powerful tool to study the influences of acquisition geometry and overburden structures on the image (e.g., Jin and Walraven, 2003; Wu and Chen, 2006; Xie et al., 2006). More accurate amplitude variation with angle (AVA) could be obtained if the image is corrected with angle dependent correction factor (e.g., Wu et al., 2004; Cao and Wu, 2005, 2008a). Traditionally, illumination analyses have used the ray-based method (e.g., Schneider and Winbow, 1999; Bear et al., 2000; Muerdter et al., 2001; Muerdter and Ratcliff, 2001a, b; Lecomte et al., 2003). The ray-based method is convenient and can provide both intensity and directional information carried in the wavefield. However, the high-frequency asymptotic approximation and the caustics inherent in ray theory may severely limit its accuracy in complex regions (e.g., Hoffmann, 2001). Furthermore, the ray method seemingly can provide angledependent illumination at any point; however this over-precise illumination map does not reflect the real behavior of the wavefield since it violates the Heisenberg uncertainty principle: The position and direction of a wavefield cannot be accurately specified simultaneously (Wu and Chen, 2006). To obtain reliable and frequency-dependent illumination, we need a wave-theory based method. One-way wave equation based propagators are widely used in migration and illumination analysis. Although they neglect multiples, they properly handle multi-forward-scattering phenomena, including focusing/defocusing, and diffraction. They can extrapolate the wavefield very fast and with accurate phase within certain angle range. This makes them suitable for seismic migration. However, the inability to provide localized-angle information prevented them from being utilized for illumination calculation. Recently developed techniques, such as LSS (e.g., Xie and Wu, 2002) and beamlet decomposition (e.g., Wu et al., 2000), can decompose the wavefield into local plane waves, which are mixed domain wavefields simultaneously localized in space and direction. These methods are independent on the extrapolator and coordinate system for extrapolation. LAD illuminations have been obtained by LSS (Xie et al., 2003, 2004, 2006), Gabor-Daubechies frame (GDF) beamlet decomposition (Wu and Chen, 2002, 2006; Cao and Wu, 2008b), and local exponential frame (LEF) beamlet decomposition (Mao and Wu, 2007; Cao and Wu, 2008b). However, the amplitude of the conventional one-way wave propagator is not accurate (e.g., Zhang et al., 2003, 2005; Wu and Cao, 2005; Cao and Wu, 2005, 2008a). The amplitude errors in the conventional one-way wave propagators have been studied by different approaches for a long time, e.g., the “true-amplitude” one-way wave equation method (Zhang, 1993; Zhang et al., 2003, 2005), multi-one-way method (Kiyashchenko et al., 2005). However, even with these corrections, the one-way propagator still cannot give accurate wavefield amplitude in complex models with sharp contrast. The numerical implementations based on the one-way wave equation with z-axis as the preferred propagation direction always have inherent limitation in wide-angle accuracy. Some one-way propagation based methods can model the 68

back-scattered waves or turning waves (e.g., Wu, 1994; Xie and Wu, 2001; Kiyashchenko et al., 2005; Jin et al. 2006; Wu and Jia, 2006; Xu and Jin, 2006; Zhang et al., 2006). However the amplitude accuracy for these propagators is still an issue due to various factors. Therefore, illumination analysis using full-wave propagators is very desirable as a complementary method. The full-wave propagator, e.g., finite-difference (FD) and finite-element method, can simulate accurate wave behavior in complex media. It provides the basis for accurate illumination analysis. Considering the expensive computation and huge storage requirement for time-domain methods and the fact that the illumination is frequency-dependent, we propose to analyze the illumination in the frequency domain. We can use frequency domain full-wave forward modeling methods (e.g., Lysmer and Drake, 1972; Marfurt, 1984; Marfurt and Shin, 1989; Jo et al., 1996; Shin and Sohn, 1998; Stekl and Pratt, 1998; Min et al., 2000) to calculate the wavefield for interested frequencies. Using 1D decomposition technique LSS, Luo et al. (2004) compared the single frequency illumination using one-way wave propagator and full-wave propagator (where single frequency full waves are extracted from the results by time-domain FD). However, artificial interference patterns exist in the illumination results from the full-wave propagator. The reason is that the decomposition was applied only along the horizontal coordinate. The downgoing and upgoing waves are mixed together, making the amplitude of the illumination incorrect and causing the interference patterns. Xie and Yang (2008) and Yang et al. (2008) proposed an illumination method which uses the full-wave FD method as the propagator and uses a time domain local slowness analysis method to determine the angle information and calculate the illumination. The method is particularly useful to provide illumination analysis for reverse time migration. Here we propose an efficient split-step method using 1D/2D decomposition to obtain LAD illumination in the frequency domain for 2D/3D full-wave propagators. Firstly we summarize the LSS and beamlet decomposition techniques for one-way propagators; based on these we give the wavefield decomposition methods for full-wave propagators. Following that we show the illumination in several models with the proposed method and compare the results with those from more expensive 2D LSS method. Finally we discuss the influence of multiples on the illumination.

7.1. Wavefield Decomposition for One-Way Propagators The local plane wave decomposition techniques, such as the LSS method and beamlet decomposition method, were applied along the horizontal coordinate(s) (e.g., Xie and Wu, 2002; Wu and Chen, 2002; Xie et al., 2003), which are appropriate to one-way propagators but not for full-wave propagators. Firstly we summarize these two decomposition methods for one-way propagators. Local slant stack for one-way propagators Figure 7.1 is a sketch showing the geometry of the LAD analysis for 2D models. In LSS, using a windowed Fourier transform along horizontal coordinate we can decompose the wavefield uz  x,   at depth z for frequency ω, obtained from any extrapolator, into local plane waves u z ( x,  ,  )   w  x  x  uz  x,   e i ( x  x )k ( x )sin  dx ,

69

(7.1)

where  is the local propagating angle with respect to the vertical direction, w is a 1D window in the horizontal direction and centered at x, k (x )   v (x ) , v ( x ) is the local velocity, and x  ( x, z ) is the location. Beamlet decomposition for one-way propagators The beamlet decomposition also provides a basis to simultaneously localize the wavefield in space and wavenumber domain. The wavefield uz  x,   can be decomposed into beamlets, uz ( x,  )   uˆz ( xn , m ,  )bmn ( x ) , m

(7.2)

n

where bmn are the decomposition basis vectors (beamlets) located at space window xn and

Figure 7.1: Basic geometry of the local angle domain analysis for 2D model (From Xie et al., 2003). wavenumber window  m , and uˆ z ( xn , m ,  ) are the decomposition coefficients for the corresponding beamlets bmn. We can obtain the local plane waves uz ( x, m ,  ) by partial reconstruction of the beamlet domain wavefield, uz ( x, m ,  )   uˆ z ( xn , m ,  )bmn ( x ) .

(7.3)

n

Both the GDF and the LEF beamlets bear uniquely defined directional localization; however, the GDF beamlet decomposition method can provide more accurate directional wavefield than the LEF beamlet decomposition (Cao and Wu, 2008b). Appendix A summarizes the GDF beamlet decomposition and partial reconstruction to obtain the local plane waves.

7.2. Wavefield Decomposition for Full-Wave Propagators When the local Fourier transform or beamlet decomposition is applied only along horizontal coordinate, the local plane wave uz ( x,  ,  ) or uz ( x, m ,  ) includes not only the wave with positive vertical wavenumber (propagating downward) but also corresponding negative vertical wavenumber (propagating upward). In one-way propagators, the waves only propagate along one primary direction; therefore decomposition techniques applied only along horizontal coordinate is appropriate. However, for the full wave propagators, they include both the down- and upgoing waves. The decomposition techniques applied only along horizontal coordinate will mix the down- and up-going waves, resulting in incorrect illumination amplitude and artificial interference patterns in the illumination map for full wave propagators (Figure 7.2, Figure 7.3; in all figures, red and blue mean strong and weak illumination respectively). 2D decomposition for full-wave propagators

70

One direct way to obtain the LAD wavefield for full waves is using 2D local plane wave decomposition for 2D problems. Similar to 1D LSS, we can decompose the wavefield for a given frequency into local plane waves using 2D LSS u ( x,  ,  )   w  x  x  u  x,   e ik ( x )(( x  x )sin   ( z  z ) cos ) dx ,

(7.4)

where w is a 2D spatial window. Split-step decomposition for full-wave propagators 2D decomposition can directly obtain the wavefield for all directions; however, it costs much more computation than 1D decomposition. Here we propose a split-step decomposition method to obtain the LAD wavefield for full-wave propagators using more efficient 1D decomposition technique. Firstly, we decompose the full-wave along vertical direction using 1D technique to separate the downgoing and upgoing waves; then apply 1D decomposition along the horizontal direction to the downgoing or upgoing waves to obtain corresponding LAD wavefield. The proposed method can obtain the illumination not only for downgoing waves, but also for the upgoing waves including turning waves and reflected waves. This method can be extended to 3D case with 1D decomposition along vertical direction and 2D decomposition along horizontal coordinates.

0

2

4

Surface location (km) 6 8 10

12

14

Depth (km)

(a) 1 2 3

Depth (km)

(b) 1 2 3

Depth (km)

(c) 1 2 3

Figure 7.2: Directional illumination for full waves with 1D LSS applied along horizontal direction for different incident angles: (a) 0º; (b) +40º; (c) -40º. Note the interference patterns (strips) caused by the interaction between up- and down-going waves.

71

0

2

4

Surface location (km) 6 8 10

12

14

Depth (km)

(a) 1 2 3

Depth (km)

(b) 1 2 3

Depth (km)

(c) 1 2 3

Figure 7.3: Acquisition dip response for full waves with 1D LSS applied along horizontal direction for different dip angles: (a) 0º; (b) +30º; (c) -30º. Note the interference patterns (strips) caused by the interaction between up- and down-going waves.

For the first step in the proposed method, we need to separate the waves with positive and negative vertical wavenumbers. Both the GDF and the LEF beamlet decomposition are very efficient to provide the local wavenumber domain wavefield with uniquely defined directional localization; hence they can be used for this step. We utilize the GDF beamlet decomposition in the first step since it can provide more accurate directional wavefield than the LEF beamlet decomposition (Cao and Wu, 2008b). For the second step, we can use either the LSS method or more efficient method with the GDF beamlet decomposition (Cao and Wu, 2008b).

7.3. Local-Angle Domain Illumination Results We exemplify the method with illumination, including DI and ADR, for downgoing waves in the 2D SEG/EAGE salt model (Aminzadeh et al., 1994; Aminzadeh et al., 1995), turning waves in a v(z) model, and reflected waves in a two-layer model. Downgoing waves illumination in 2D SEG/EAGE salt model The whole acquisition system of this model consists of 325 shots with 176 left-side trailing receivers for each shot and the shot and receiver intervals are 160 and 80 feet respectively. The dominant frequency for the source wavelet is 15 Hz. In the following, we will only consider the dominant frequency for both DI and ADR calculations. DI results for incident angles 0º, +40º and -40º from the proposed method are similar to those obtained by the much more expensive 2D LSS method (Figure 7.4). DI for incident angles from 60º to +60º at different location along depth level z=3.2675 km shows that these LAD wavefield 72

obtained by the proposed split-step decomposition method are indistinguishable from those by the 2D LSS method (Figure 7.5). For ADR, the results for dip angles 0º, +30º and -30º from the proposed method are also similar to those obtained by the 2D LSS method (Figure 7.6). ADR for dip angles from -60º to +60º at different location along depth level z=3.2675 km also show that the result by the split-step method is almost the same as that by the 2D LSS method (Figure 7.7).

0

2

4

Surface location (km) 6 8 10

12

14

Depth (km)

(a) 1 2 3

Depth (km)

(b) 1 2 3

Depth (km)

(c) 1 2 3

Depth (km)

(d) 1 2 3

Depth (km)

(e) 1 2 3

Depth (km)

(f) 1 2 3

Figure 7.4: Directional illumination for full waves from the proposed method (a-c) and 2D LSS (d-f) for different incident angles: (a, d) 0º; (b, e) +40º; (c, f) -40º.

73

Incident angle (deg)

-60 -30 0 30 60 -60 -30 0 30 60

0

2

4

Surface location (km) 6 8 10

12

14

(a)

(b)

Figure 7.5: Directional illumination for incident angles from -60º to +60º at different location along depth level z=3.2675 km: (a) from the proposed method; (b) from 2D LSS. 0

2

4

Surface location (km) 6 8 10

12

14

Depth (km)

(a) 1 2 3

Depth (km)

(b) 1 2 3

Depth (km)

(c) 1 2 3

Depth (km)

(d) 1 2 3

Depth (km)

(e) 1 2 3

Depth (km)

(f) 1 2 3

Figure 7.6: Acquisition dip response for full waves from the proposed method (a-c) and 2D LSS (d-f) for different dip angles: (a, d) 0º; (b, e) +30º; (c, f) -30º. 74

Dip angle (deg)

-60 -30 0 30 60 -60 -30 0 30 60

0

2

Surface location (km) 6 8 10

4

12

14

(a)

(b)

Figure 7.7: Acquisition dip response for dip angles from -60º to +60º at different location along depth level z=3.2675 km: (a) from the proposed method; (b) from 2D LSS. Turning waves and reflected waves illumination Turning waves and reflected waves can well image the overhung or vertical structures (e.g., salt flanks), which the downgoing waves in traditional one-way propagators cannot image (e.g., Jin et al. 2006; Xu and Jin, 2006; Zhang et al., 2006; Jia and Wu, 2007). Here we illustrate the illumination for the turning waves and reflected waves obtained from the full-wave propagator. The examples are for one shot and frequency 15Hz. Figure 7.8 shows the DI for incident angle 135º in a v(z) model. The result from the proposed method is very similar to that obtained by 2D LSS. We can also notice that the result from the proposed method has a lower resolution than that from the 2D LSS method. The artificial interference pattern is caused by the interaction of turning waves and a weak upgoing reflected wave. This reflection is produced by the velocity gradient change at the model bottom where we pad the v(z) model with a constant velocity in the full-wave modeling. We use a simple model to generate the reflected waves. It consists of a 5 km thick homogeneous layer with velocity 3.0 km/s and a half space with velocity 4.5 km/s. Figure 7.9 shows the DI for incident angle 135º. The result from the proposed method is very similar to that obtained by 2D LSS. The artificial interference pattern in the lower-right corner of the model is caused by the interaction of reflected waves and head waves.

Depth (km)

Depth (km)

0

2

4

6

Surface location (km) 8 10 12

14

16

18

1 2

(a)

3 4

1 2

(b)

3 4

Figure 7.8: Directional illumination for incident angle 135º in a v(z) model: (a) from the proposed method; (b) from 2D LSS. The star represents the source location.

75

Depth (km)

Depth (km)

0

2

4

6

Surface location (km) 8 10 12

14

16

18

1 2

(a)

3 4

1 2

(b)

3 4

Figure 7.9: Directional illumination for incident angle 135º in a two-layer model: (a) from the proposed method; (b) from 2D LSS. The star represents the source location.

7.4. Discussion About The Influence of Multiples on The Illumination Strength Above proposed split-step decomposition method can separate the downgoing and upgoing waves. The downgoing waves include not only the primary incident waves but also multiples; and the primary incident waves include the first arrival and multi-arrivals. Figure 7.10 shows the DI maps for the most energetic waves. Comparison with the DIs for the full downgoing waves (Figure 7.4a-c) shows that the other arrivals (multi-arrivals and multiples) provide extra illumination to the subsurface. The migration methods based on one-way propagators utilize not only the first arrival but also the multi-arrivals. Therefore the multi-arrivals should be included in the illumination analysis for survey design and true-reflection imaging correction. However, multiples especially the internal multiples should be eliminated since most of the migration methods based on the one-way propagator only utilize primary incident waves. It is impossible to separate multiples from other arrivals in the frequency domain. In time domain, we can separate the first arrival or the most-energetic arrival from other arrivals by time-domain windowing. However, we cannot separate multiples from multi-arrivals by time windowing since they may have similar travel time. In general media, it is hard to obtain the pure multiple data or remove them from the full data. The one-way and one-return boundary element method in the frequency domain (He and Wu, 2007) can calculate the primary transmitted waves and multiples for layered model or inclusion model. This method can handle strong velocity contrast. Here we exemplify the influence of internal multiples on the illumination strength with a simplified SEG/EAGE salt model. A homogeneous salt body with the same shape as that in the original 2D SEG/EAGE salt model is embedded in a homogeneous background medium. The velocity of the background and the salt are

76

0

2

Surface location (km) 6 8 10

4

12

14

Depth (km)

(a) 1 2 3

Depth (km)

(b) 1 2 3

Depth (km)

(c) 1 2 3

Figure 7.10: Directional illumination for the most energetic waves with the proposed method for different incident angles: (a) 0º; (b) +40º; (c) -40º.

Depth (km)

(a)

6

(b)

12

1.8

1.8

2.6

2.6

3.4

6

Surface location (km) 8 10

12

3.4 0

(c) Depth (km)

Surface location (km) 8 10

0.2 6

0.4 8

0

0.6 10

(d)

12

1.8

1.8

2.6

2.6

3.4

0.02 6

0.04 8

10

0.06 12

3.4 0

0.1

0.2

0.3

0

0.4

0.02

0.04

0.06

Figure 7.11: Directional illumination for the transmitted waves and salt internal multiples in the simplified SEG/EAGE salt model: (a) transmitted waves for -10º; (b) multiples for -10º; (c) transmitted waves for -60º; (d) multiples for -60º. 7800 ft/s and 14700 ft/s respectively. This is a scalar wave model. For acoustic wave model, the salt internal multiples will be weaker since the density of the salt is usually lower than that of the background sediment. DI maps in sub-salt region for incident angles -10º and -60º for the primary incident waves and salt internal multiples (Figure 7.11) show that the maximum amplitude of the primaries is more than six times stronger than that of the multiples for these two angles. Therefore the contribution of the multiples to the illumination strength could be considered as a secondary effect here. However, we can also notice the difference in the spatial distribution of the illumination between the results for multiples and those for primaries. For example, for -60º incident angle, the primaries only strongly illuminate the left part of the sub-salt; however, the 77

multiples illuminate the whole sub-salt more evenly and provide extra illumination to the shadow in the illumination by primaries (right part in sub-salt).

7.5. Conclusions We propose an efficient split-step local plane wave decomposition method to obtain the localangle domain (LAD) illumination in the frequency domain for full-wave propagators. It can obtain the illumination not only for the downgoing waves but also for the upgoing waves including turning waves, refracted and reflected waves. It can be used for accurate survey design and true-reflection imaging correction to obtain accurate AVA for media property inversion. Numerical examples show that directional illumination and acquisition dip response results from the proposed method are very similar to those obtained by the more expensive 2D decomposition method, although the proposed method could has a lower resolution. The results in a simplified SEG/EAGE salt model show that: the illumination strength from salt internal multiples could be considered as a secondary effect compared with that from primary incident waves; however, the multiples can illuminate the whole sub-salt more evenly and provide extra illumination to the shadow in the illumination by primaries. The proposed method has following features: 1. It is in the frequency domain, so it can provide frequency-dependent illumination and the frequency-domain full-wave forward modeling is usually efficient and storage-saving to provide the wavefield for given frequencies compared with the time-domain method; 2. It is efficient compared with the full 2D/3D decomposition since we only need three 1D decompositions for 2D models and two 2D plus one 1D decompositions for 3D models; 3. It is independent of the propagator used for full-wave modeling. We conclude that the proposed split-step local plane wave decomposition method for full-wave propagators provides a flexible, efficient and accurate tool for the LAD wave-theory based illumination analysis in complex 2D and 3D models.

Appendix A Gabor-Daubechies Frame (GDF) Beanlet Decomposition The GDF beamlets (e.g., Wu et al., 2000; Wu and Chen, 2001; Chen et al., 2006) have uniquely defined and very good localization information available after the decomposition. The GDF beamlets for the beamlet decomposition in equation (7.2) are Gaussian function windowed exponential harmonics, bmn ( x )  g ( x  xn )eim x , (A-7.1)

where m  m (  is the wavenumber sampling interval), and g  x  is a Gaussian window function,  x2  (A-7.2) g ( x )  ( s 2 ) 1/ 4 exp   2  ,  2s  s2 

R  N 2 , 2

(A-7.3)

78

where s is the scale of the Gaussian window, R is the redundancy ratio and N is the lateral sampling interval of the frame. Substituting equation (A-7.1) into (7.3), we can obtain the partially reconstructed local plane waves, (A-7.4) u z ( x,  m ,  )  eim x  g ( x  xn )uˆ z ( xn ,  m ,  ) , n

with the decomposition coefficients * uˆ z ( xn , m ,  )  uz ( x,  ), bmn ( x )   uz ( x,  )  bmn ( x ) dx ,

(A-7.5)

where , stands for inner product, “*” stands for complex conjugate, and bmn ( x ) are the dual GD frame atoms, bmn ( x )  g ( x  xn )eim x , (A-7.6) with g ( x ) as the dual-window function of g ( x ) . The dual-window function can be calculated by pseudo-inversion of the original window function (Qian and Chen, 1996; Mallat, 1998; Wu and Chen, 2001). From equation (A-7.4), we can see that the local plane wave is a weighted average of the windowed beamlets with same wavenumber from neighboring windows. The total space domain wavefield can be written as (A-7.7) uz ( x,  )   uz ( x,  m ,  ) . m

79

8. FAST ACQUISITION APERTURE CORRECTION BY BEAMLET DECOMPOSITION Seismic imaging should provide us an estimate of the subsurface reflectivity/scattering strength, which is reflection/scattering angle dependent. Therefore the imaging and image amplitude correction have to be done in the angle-domain, which makes ray-based methods more convenient than wave-equation based methods since the angle information is inherently embedded in ray-based methods. The theory and method of true-reflection imaging has been developed based on high-frequency asymptotic theory (ray theory) and is traditionally carried out through Kirchhoff prestack depth migration (e.g., Bleistein et al., 1987; Hubral et al., 1991; Hanitzsch, 1995; Xu et al., 2001; Audebert et al., 2002; Brandsberg-Dahl et al., 2003). However, the results may contain large errors in complex areas due to the high-frequency approximation and singularity problems therein. For wave-equation based migration methods, Wu et al. (2004) proposed an image amplitude correction method in LAD for acquisition aperture. Numerical examples showed significant improvement in both the total strength of the images and angle-dependent reflection amplitudes. This demonstrated the significance of acquisition aperture correction in true-reflection imaging. For wave-equation based migration methods, the subsurface angle-domain common-image gathers (CIGs) can be obtained either by decomposing the propagated wavefield before imaging (“data-space methods”) (e.g., de Bruin et al., 1990; Mosher et al., 1997; Prucha et al., 1999; Mosher and Foster, 2000; Wu and Chen, 2002; Xie and Wu, 2002; Chen et al., 2006), or by prestack imaging at zero time but not at zero offset (“image-space methods”) (e.g., Rickett and Sava, 2001, 2002; Biondi and Shan, 2002; Sava and Fomel, 2003; Biondi and Symes, 2004; Biondi 2007; Rosales et al, 2008). In the image-space method, the subsurface offset-domain CIGs are firstly computed and then transformed to the subsurface angle-domain CIGs. For source-receiver migration, the subsurface offset-domain CIGs are naturally available. For shot profile migration, the subsurface offset-domain CIGs are obtained by crosscorrelating the source and receiver wavefields shifted horizontally with respect to each other extended from conventional crosscorrelation imaging condition in time (Rickett and Sava, 2001, 2002). The prestack image becomes a function of the horizontal relative shift, which has the physical meaning of subsurface half offset. As the geologic dips increase, the subsurface horizontal-offset CIGs degenerate, and their focusing around zero offset blurs (Biondi and Symes, 2004). This limitation of subsurface horizontal-offset CIGs can be sidestepped by computing subsurface offset-domain CIGs along the vertical subsurface-offset axis (Biondi and Shan, 2002; Biondi and Symes, 2004). Sava and Fomel (2003) presented a simple method for transforming subsurface horizontal-offset CIGs into subsurface angle-domain CIGs by a slant stack transformation (Schultz and Claerbout, 1978) applied independently to each subsurface horizontal-offset CIG. Biondi and Symes (2004) proposed a simple and effective method to create a subsurface angle-domain CIG cube by combining a subsurface horizontal-offset CIG cube with a subsurface vertical-offset CIG cube. In data-space methods, angle gathers are obtained by decomposing the propagated wavefield before imaging. The angle-domain images can be obtained for shot-profile migration by de Bruin et al. (1990) using global Fourier transform in 1D media, and for source-receiver migration (Mosher et al., 1997; Prucha et al., 1999; Mosher and Foster, 2000). Recently developed 80

techniques, such as the local slant stack method proposed by Xie and Wu (2002) and the beamlet decomposition method proposed by Wu and Chen (2002), can decompose the wavefield into local plane waves, which are simultaneously localized in space and direction. These methods are more general and have advantages: 1. they are independent of the extrapolator (no matter it is shot-profile, source-receiver, or other propagators); 2. with the decomposed local plane waves we can obtain not only the angle-domain image but also many other very useful angle-domain information, e.g., LAD illumination (Wu and Chen, 2002, 2006; Xie et al., 2003, 2004, 2006), resolution (Xie et al, 2005; Wu et al, 2006), and image amplitude compensation factor (Wu et al., 2004). The implementation of the original acquisition aperture correction (Wu et al., 2004) is very timeconsuming since the decomposition of the wavefield uses LSS which is computationally demanding. For example, even using a relatively short spatial window for decomposition, the computation time to obtain the image in local dip-angle domain can be five times of that to obtain the space-domain image. This prevents the acquisition aperture correction in LAD from practical application in industrial migration processing. Beamlet decomposition (e.g., Wu et al., 2000) can also decompose the wavefield into local plane waves. The propagation and imaging methods based on beamlet decomposition proposed by Wu et al. (2000) have been developed using the GDF (Wu et al., 2000; Wu and Chen, 2001; Chen et al., 2006), LCB (local cosine basis) (Wu et al., 2000; Wang and Wu, 2002; Luo and Wu, 2003; Wang et al., 2003; Wu et al., 2003; Wu et al., 2008), and LHB (local harmonic basis) (Wu and Mao, 2007). The GDF beamlets have uniquely defined local direction information available during the propagation process so that angle-domain imaging (Wu and Chen, 2002; 2003; Chen et al., 2006) and other angle-related operations can be easily incorporated into the migration process, but the GDF beamlet propagator is computationally more expensive than the LCB propagator due to the redundancy in its representation. Although LCB beamlets are orthonormal, they always have two symmetrical lobes with respect to the vertical axis due to the inherent property of the cosine basis function. This lack of uniquely defined directional localization prevents its use for operations in the LAD. Local exponential frame (LEF) (Daubechies et al., 1991; Auscher, 1994; Mao and Wu 2007) beamlet decomposition can eliminate the directional symmetry in the LCB beamlets. It can be implemented by a combination of local cosine and sine transforms which have fast algorithms; therefore the LEF decomposition is efficient. Beamlet decomposition is more efficient than LSS; however the obtained wavefield is in the LWD not directly in the LAD. Original LAD imaging with beamlet decomposition transforms the complex wavefield from LWD to LAD, which is also very time consuming. We propose a fast implementation to obtain the LAD image and amplitude correction factor by beamlet decomposition in the LWD. Firstly the beamlet decomposition method is summarized. Then the new implementation is described. We exemplify the method with the GDF and the LEF decomposition in 2D SEG/EAGE salt model (Aminzadeh et al., 1994; Aminzadeh et al., 1995) and compare the results with those from the original method. We also discuss the possible factors causing the difference in the results by the GDF and the LEF decomposition.

8.1. Beamlet Decomposition of the Wavefield

81

The beamlet decomposition (e.g., Steinberg, 1993; Wu et al., 2000; Wu and Chen, 2001; Chen et al., 2006; Mao and Wu, 2007) provides a basis to simultaneously localize the wavefield in space and wavenumber domain. The frequency-space domain wavefield u ( x,  ) at depth level z for frequency  , obtained from any extrapolator, can be decomposed into beamlets, u( x,  )   uˆ ( xn ,  m ,  )bmn ( x ) , (8.1) m

n

where bmn are the decomposition basis vectors (beamlets) located at space window xn and wavenumber window  m , and uˆ( xn ,  m ,  ) are the decomposition coefficients for corresponding beamlets bmn. We can obtain the local plane waves u( x, m ,  ) (mixed domain wavefield in local phase space) by partial reconstruction of the beamlet domain wavefield, u( x,  m ,  )   uˆ ( xn ,  m ,  )bmn ( x ) . (8.2) n

The beamlet used to obtain the LAD image and image amplitude correction factor must have uniquely defined local direction information. In this paper, we exemplify the new method with two types of beamlet decomposition: the GDF beamlet and the LEF beamlet. Appendix A and B summarize the GDF and LEF beamlet decomposition and partial reconstruction to obtain the local plane waves respectively. Both the GDF and the LEF beamlet decomposition methods use windowed exponential harmonics as the decomposition basis. The differences are the windows used to modulate the harmonics and the redundancy ratio (ratio of the number of all decomposed beamlet coefficients to the length of the signal in space domain, which describes the overlapping among neighboring frames). Bell window (see, e.g., Wu et al., 2008) and Gaussian window (see, e.g., Wu and Chen, 2006) are used in the LEF and GDF decomposition respectively. The redundancy ratio R in the GDF decomposition is adjustable, which make it be able to provide very accurate partially reconstructed LWD wavefield. However, R in the LEF representation is fixed to 2. It could cause some error in the LWD wavefield, which will be shown in the paper. For local plane wave with local wavenumber  m , the corresponding propagating angle with respect to the vertical direction is  m  sin 1 m k (x )  , (8.3)

where k ( x )   v ( x ) and v ( x ) is the local velocity.

8.2. LAD Image and Image Amplitude Correction Factor by LWD Bemlet Decomposition For subsurface imaging point, the traditional space-domain imaging condition (Claerbout, 1971) for a single frequency is in the form of crosscorrelation of the incident source-side wave uS ( x,  ) and the back-propagated receiver-side recorded wave uRs (x,  ) , I ( x,  )   uS* ( x,  )  uRs ( x,  ) S

  G ( x; x s )  2  * I

,

(8.4)

d (x g ; x s ) z where “*” stands for complex conjugate; GI is the Green’s function used in the imaging process; and the integral is a back propagation Rayleigh integral of the recorded data, in which A(xg; xs) is xs

A ( x g ;x s )

82

dx g

GI* ( x; x g )

the spatial receiver aperture for given source, d (xg; xs) is the recorded scattered data at the receiver xg from the source at xs. Seismic imaging should provide us an estimate of the subsurface reflectivity/scattering strength, which is reflection/scattering angle dependent. Therefore the conventional imaging condition in the space domain (8.4) needs to be extended to the space-local angle domain (or beamlet domain) (e.g., Wu and Chen, 2002, 2003; Chen et al., 2006). Then the image obtained at each imaging point is no longer a scalar but a matrix, called the local image matrix (LIM) L(x, θs, θg), where θs and θg are the source and receiving angles, respectively (Figure 8.1). The LIM is a distorted estimate of the local scattering matrix (LSM) due to the acquisition aperture limitation and the propagation path effects. LSM is the intrinsic property of the scattering medium and is independent of the acquisition system and free from propagation effects and contains information of the local structure and elastic properties (Wu et al., 2004). The task of true-reflection imaging is to restore the true LSM from the distorted LIM by applying image amplitude corrections. In the original LAD imaging implementation, the image is directly formed using the wavefield in the LAD, which can be obtained by either LSS (e.g., Xie and Wu, 2002; Wu et al., 2004) or beamlet decomposition plus complex interpolation. Both methods are very time-consuming. In the original LAD implementation with beamlet decomposition, for the local plane wave with

Figure 8.1: Local-angle definition in a scattering experiment for a local planar reflector. θs and θg are the local source and receiving angles, respectively. θn is reflector-normal angle (equal to dipangle) and θr is the reflection angle. direction other than the discretely sampled lobe direction  m , interpolation between neighboring lobes is used. This interpolation is applied to the complex incident and back-propagated wavefields for every shot gather; therefore, it is very time-consuming. For LSS method, the decomposition even has to be applied at every imaging point. Since it is much more efficient to obtain the wavefield in the LWD than in the LAD, we propose an efficient way to obtain the LAD image by partial imaging in the LWD: firstly the image in the LWD for a single frequency can be formed for each shot and then summed up to form the LWD total image for all shots which is then transformed to the LAD. Using the local plane waves in the LWD obtained by beamlet decomposition (e.g., equation (8.2)), the imaging condition for a single frequency ω in the space-local wavenumber domain can be written as, L( x,  m ,  j ,  )   uS* ( x,  m ,  )  uRs ( x,  j ,  ) , (8.5) S

83

where  m and  j are the source and receiving wavenumbers respectively; uS ( x,  m ,  ) is the incident beamlet at the imaging point x generated by the source S at xs; and uRs (x,  j ,  ) is the back-propagated beamlet at the imaging point x from the recorded data for source S. The summed LWD image for all shots and a given frequency L(x,  m ,  j ,  ) in (8.5) is a real number and can easily be interpolated to form the final LAD image L(x, θs, θg, ω). Since we only need to transform this shot-summed real image matrix from the LWD to the LAD by interpolation, the computation time to obtain the final LAD image for all shots is significantly reduced. The LAD images from different frequencies can then be summed to form the multi-frequency image. With similar derivation to that obtaining the image amplitude correction factor directly from the LAD wavefield (Wu et al., 2004), we can obtain the LWD image amplitude correction factor for corresponding imaging condition (8.5), Fa (x,  m ,  j )  2 G ( x, m ; x s )GF ( x,  m ; x s )  * I

xs



A( x g ;x s )

dx g GF ( x,  j ; x g )

2



1/2

,

(8.6)

where GF is the Green’s function used in forward modeling. With similar scheme to obtaining the LAD image in previous part, the LWD amplitude correction factor (8.6) after shot summation can also be interpolated to that in the LAD, Fa(x, θs, θg). Having obtained the LAD image matrix L(x, θs, θg) and amplitude correction factor matrix Fa(x, θs, θg), we can apply same types of acquisition aperture correction as those in the original method, such as the correction for common reflection-angle imaging and correction for total strength imaging.

8.3. Examples We exemplify the new implementation of LAD acquisition aperture correction with 2D SEG/EAGE salt model (Figure 8.2) data set (Aminzadeh et al., 1994; Aminzadeh et al., 1995). The whole acquisition system of this model consists of 325 shots with 176 left-side trailing receivers for each shot and the shot and receiver intervals are 160 and 80 feet respectively. We employ LCB propagator (Wu et al., 2000; Wang and Wu, 2002; Luo and Wu, 2003; Wang et al., 2003; Wu et al., 2003, 2008) for the wavefield extrapolation. The redundancy ratio R in the GDF decomposition is set to 4 except specified in the following numerical examples. X

The images in the local dip-angle domain obtained from the new implementation by the LEF decomposition (Figure 8.3) and the GDF decomposition (Figure 8.4) are very similar to those from the original method with LSS (Figure 8.5). Some subtle differences exist among these images from different methods, e.g., the artifacts around the steep salt top in the images obtained from the LEF decomposition are stronger than those from the other two methods, and the images from the GDF decomposition are almost the same as those from the original method. The image amplitude correction factors in the local dip-angle domain from the new method by the LEF decomposition (Figure 8.6) ; in the figures, red and blue mean strong and weak factor amplitude respectively) and the GDF decomposition (Figure 8.7) are also very similar to those from the original method (Figure 8.8). The results from the new method by the GDF decomposition are almost the same as the original ones; however, in the results obtained by the LEF decomposition X

84

there are some stripe-like artifacts (the width of the stripe is related to the nominal length of the window (or frame sampling interval) in the LEF decomposition). We will discuss the reasons causing these artifacts in following “Discussion” part. The final total strength of the image after acquisition aperture correction in the local dip-angle domain from both the LEF and GDF decomposition is also very similar to that from the original method (Figure 8.9). Similar to the images before correction, there are also some subtle differences among the final corrected images: the artifacts in the subsalt area and around the steep salt top obtained from the LEF decomposition are stronger than that from the other two methods, and the corrected image by the GDF decomposition is almost the same as the original one. As for the efficiency, if we set the computation cost of the conventional space-domain imaging to 1, the cost of the original method to form the LAD image with LSS is about 5.2 even using a relatively short spatial window for decomposition. For the new method it is only about 1.9 for the LEF decomposition and 2.1 for the GDF decomposition, which is more than twice faster than the original method and only costs about twice computation time of the traditional space-domain imaging. Here we employ LCB beamlet propagator for the extrapolation, but the computation time by the LEF decomposition is obtained without taking advantage of the availability of local cosine coefficients during LCB beamlet propagation. If taking that advantage it will save half of the time for obtaining the beamlet decomposition coefficients although it is only a very small fraction in the total computation time.

Depth (x 80 feet)

0

0

Surface location (x 80 feet) 200 300 400

100

500

600

50 100

Figure 8.2: 2D SEG/EAGE salt model.

0

0

100

200

300

400

500

600

50

(a) 100

Image in dip angle domain (0 degree)

85

0

0

100

200

300

400

500

600

50

(b) 100

Image in dip angle domain (+30 degree) 0

0

100

200

300

400

500

600

50

(c) 100

Image in dip angle domain (-30 degree)

Figure 8.3: The image obtained from the new method by the LEF decomposition for different dip angles: (a) 0º; (b) +30º; (c) -30º.

86

0

0

100

200

300

400

500

600

50

(a) 100

Image in dip angle domain (0 degree)

0

0

100

200

300

400

500

600

50

(b) 100

Image in dip angle domain (+30 degree)

0

0

100

200

300

400

500

600

50

(c) 100

Image in dip angle domain (-30 degree)

Figure 8.4: The image obtained from the new method by the GDF decomposition for different dip angles: (a) 0º; (b) +30º; (c) -30º.

87

0

0

100

200

300

400

500

600

50

(a) 100

Image in dip angle domain (0 degree) 0

0

100

200

300

400

500

600

50

(b) 100

Image in dip angle domain (+30 degree)

0

0

100

200

300

400

500

600

50

(c) 100

Image in dip angle domain (-30 degree)

Figure 8.5: The image obtained from the original method by LSS for different dip angles: (a) 0º; (b) +30º; (c) -30º.

88

0

0

100

200

300

400

500

600

50

(a) 100

Factor in dip angle domain (0 degree) -new method 0

0

100

200

300

400

500

600

50

(b) 100

Factor in dip angle domain (+30 degree) -new method

0

0

100

200

300

400

500

600

50

(c) 100

Factor in dip angle domain (-30 degree) -new method

Figure 8.6: The image amplitude correction factor obtained from the new method by the LEF decomposition for different dip angles: (a) 0º; (b) +30º; (c) -30º.

89

0

0

100

200

300

400

500

600

50

(a) 100

Factor in dip angle domain (0 degree) -new method (GDF)

0

0

100

200

300

400

500

600

50

(b) 100

Factor in dip angle domain (+30 degree) -new method (GDF)

0

0

100

200

300

400

500

600

50

(c) 100

Factor in dip angle domain (-30 degree) -new method (GDF)

Figure 8.7: The image amplitude correction factor obtained from the new method by the GDF decomposition for different dip angles: (a) 0º; (b) +30º; (c) -30º.

90

0

0

100

200

300

400

500

600

50

(a) 100

Factor in dip angle domain (0 degree) -original method 0

0

100

200

300

400

500

600

50

(b) 100

Factor in dip angle domain (+30 degree) -original method

0

0

100

200

300

400

500

600

50

(c) 100

Factor in dip angle domain (-30 degree) -original method

Figure 8.8: The image amplitude correction factor obtained from the original method by LSS for different dip angles: (a) 0º; (b) +30º; (c) -30º.

91

0

0

100

200

300

400

500

600

50 100

(a)

Before correction

0

0

100

200

300

400

500

600

50

(b) 100

Fast acquisition aperture correction (LEF)

0

0

100

200

300

400

500

600

50

(c) 100

Fast acquisition aperture correction (GDF)

0

0

100

200

300

400

500

600

50

(d) 100

Original acquisition aperture correction (LSS)

Figure 8.9: The final total strength of the image after acquisition aperture correction in the local dip-angle domain: (a) without correction; (b) correction with the new fast method by the LEF decomposition; (c) correction with the new fast method by the GDF decomposition; (d) correction with the original method by LSS.

92

8.4. Discussion Above examples show that the results from the GDF decomposition are superior to those from the LEF decomposition, e.g., there are no stripe-like artifacts which appear in the amplitude correction factor obtained by the LEF decomposition. The differences in the LEF and the GDF beamlet decomposition are the window functions (bell window vs. Gaussian window, see (Figure 8.10) and redundancy ratio. The redundancy ratio is fixed at 2 in the LEF decomposition, but is adjustable in the GDF decomposition (4 used in our calculations). Here we investigate these two factors to find what may cause the difference in above results. We use the directional illumination (see, e.g. Wu and Chen, 2002; Xie et al., 2003) in the LWD from a single shot (the fifth shot here) for illustration because it represents the amplitude of the decomposed directional wavefield. We also investigate the amplitude correction factor for the whole acquisition system. First, we examine the influence of the window function. In previous calculations with the LEF decomposition, a relatively steep bell function with frame sampling interval N  8 (solid line in Figure 8.10) is used. With a less steep bell function (dotted line in Figure 8.10), the single shot directional illumination for the incident wave close to vertical direction shows some improvement compared with the result using the original window (Figure 8.11). There is no vertical beamlet in the LEF decomposition (see, e.g., Appendix B), so we show the results for the beamlet closest to the vertical direction. Results also show that the DI gets smoother with larger N (Figure 8.11). With the new window, the stripe-like artifacts are removed in the amplitude correction factor for dip-angle +30º in the subsalt and right part of the model (Figure 8.12a); however, those in the left part of the model are still similar to those with the steep bell function in Figure 8.6. With a larger N (=16), the artifacts are weaker, but we also lose the spatial resolution in some degree in the subsalt area (Figure 8.12b). These results demonstrate that with a longer and less steep bell window the artifacts in the LAD results by the LEF decomposition get weaker but cannot be completely removed. Next we investigate the effect of the redundancy ratio in the GDF decomposition. X

X

X

X

X

Fi

X

Amplitude

1

0.5

0

2

4

6

8

10

12

14

16

Figure 8.10. Windows used in the LEF and the GDF beamlet decomposition. Solid line: steep bell window function used in previous LEF decomposition; dotted line: smooth bell window function for the LEF decomposition; dashed line: Gaussian window function used in the GDF decomposition.

93

0

0

100

200

0

100

200

50

(a)

(b)

(c)

(d)

100 0 50 100

Figure 8.11: Directional illumination maps for the incident wave closest to the vertical direction from the 50th shot in the SEG/EAGE salt model by the LEF decomposition using bell windows with different steepness and lateral sampling interval of the frame ΔN: (a) with a less steep bell window (dotted line in Figure 8.10) and ΔN=8; (b) with same less steep bell window but ΔN=16; (c) with a steep bell window (solid line in Figure 8.10) and ΔN=8; (d) with same steep bell window but ΔN=16. X

X

0

0

100

200

300

400

500

600

50

(a) 100

Factor in dip angle domain (+30 degree) -new method 0

0

100

200

300

400

500

600

50

(b) 100

Factor in dip angle domain (+30 degree) -new method

Figure 8.12: The image amplitude correction factor for dip-angle +30º obtained from the new method by the LEF decomposition with a less steep bell window function (dotted line in Figure 8.10) but different lateral sampling interval of the frame ΔN: (a) ΔN=8; (b) ΔN=16. X

In previous calculations with the GDF decomposition, the redundancy ratio R used is 4. With R=4, the single shot directional illumination maps for the vertically incident wave are very clean

94

(Figure 8.13a, b). With a smaller redundancy ratio R=2, the illumination results (Figure 8.13c, d) show very subtle stripe-like artifacts which are much weaker than those by the LEF decomposition (Figure 8.11). With same frame sampling interval N  4 as before (e.g., in Figure 8.7) but a smaller redundancy ratio R=2, the amplitude correction factor for dip-angle +30º by the GDF decomposition (Figure 8.14a) now displays similar stripe-like artifacts in the left part of the model to those by the LEF decomposition (Figure 8.6b, Figure 8.12). However, these artifacts are much weaker for the GDF decomposition. The artificial stripes are narrower than that by the LEF decomposition (Figure 8.6b) because N for the LEF decomposition is 8 there. Using N  8 in the GDF decomposition, the result (Figure 8.14b) shows artifacts with same width as those by the LEF decomposition using either steep (Figure 8.6b) or smooth window (Figure 8.12a), but with much weaker strength. These results demonstrate that: the GDF decomposition with lower redundancy ratio (e.g., R=2) could give LAD illumination and amplitude correction factor with very subtle stripe-like artifacts; however, with higher redundancy ratio (e.g., R=4) it could yield very neat results. X

Similar to the results by the LEF decomposition, we can also notice that the amplitude correction factor by the GDF decomposition in Figure 8.14 and Figure 8.7b has different spatial resolution in subsalt area. Single shot directional illumination also shows this feature (Figure 8.13). It is caused by the use of different parameter R and N . Different combination of R and N may give differently scaled Gaussian windows (see equation (A-2) and (A-3)) and the scale determines the beamlet resolution in the space and wavenumber/angle domain. This can also explain the resolution difference in subsalt area among the results from the LEF (Figure 8.6), the GDF (Figure 8.7), and LSS (Figure 8.8). X

0

0

100

200

X

0

X

100

X

200

50

(a)

(b)

(c)

(d)

100 0 50 100

Figure 8.13: Directional illumination maps for the vertically incident wave from the 50th shot in the SEG/EAGE salt model by the GDF decomposition with different redundancy ratio R and lateral sampling interval of the frame ΔN: (a) R=4, ΔN=4; (b) R=4, ΔN=8; (c) R=2, ΔN=4; (d) R=2, ΔN=8.

95

0

0

100

200

300

400

500

600

50

(a) 100

Factor in dip angle domain (+30 degree) -new method (GDF) 0

0

100

200

300

400

500

600

50

(b) 100

Factor in dip angle domain (+30 degree) -new method (GDF)

Figure 8.14: The image amplitude correction factor for dip-angle +30º obtained from the new method by the GDF decomposition with redundancy ratio R=2 but different lateral sampling interval of the frame ΔN: (a) ΔN=4; (b) ΔN=8.

8.5. Conclusions We propose a fast method by shot summation in the local wavenumber domain (LWD) to obtain the local-angle domain (LAD) image and image amplitude correction factor with efficient LWD beamlet decomposition. With obtained LAD image and amplitude correction factor from the new method, same types of acquisition aperture correction can be applied as in the original method. The results in the 2D SEG/EAGE salt model illustrate that the new method with both the local exponential frame (LEF) and Gabor-Daubechies frame (GDF) beamlet decomposition produces similar results to the original method. As for the efficiency, the new method is more than twice faster than the original method and only costs about twice computation time of the traditional space-domain imaging. This makes it possible for the wave-equation based space local-angle domain processing, e.g., illumination analysis and acquisition aperture correction, to be applied to the huge amount of data in industrial processing. We notice that the results from the LEF decomposition are inferior to those from the GDF decomposition although it could be a little more efficient, and the results from the GDF decomposition can be artifacts free and indistinguishable from those obtained by the original method. We investigate two possible reasons, decomposition window and redundancy ratio, causing the difference in the results. Our analysis demonstrates that: with lower redundancy ratio (R=2), the GDF decomposition with Gaussian window could give LAD illumination and amplitude correction factor with stripe-like artifacts similar to those by the LEF decomposition but with much weaker strength; however it could yield very clean results with higher redundancy ratio (e.g., R=4); for the LEF decomposition, even with a longer and less steep bell window the artifacts in the LAD results

96

cannot be completely removed although they get weaker. This may be attributed to its fixed redundancy ratio 2.

Appendix A: Gabor-Daubechies Frame (GDF) Beanlet Decomposition The GDF beamlets (e.g., Wu et al., 2000; Wu and Chen, 2001; Chen et al., 2006) have uniquely defined and very good localization information available after the decomposition. The GDF beamlets for the beamlet decomposition in equation (8.1) are Gaussian function windowed exponential harmonics, bmn ( x )  g ( x  xn )eim x , (A-8.1)

where  m  m (  is the wavenumber sampling interval), and g  x  is a Gaussian window function,

 x2  (A-8.2) g ( x )  ( s 2 ) 1/4 exp   2  ,  2s  R  N 2 s2  , (A-8.3) 2 where s is the scale of the Gaussian window, R is the redundancy ratio and N is the lateral sampling interval of the frame. Substituting equation (A-1) into (8.2), we can obtain the partially reconstructed local plane waves, u( x,  m ,  )  eim x  g ( x  xn )uˆ ( xn ,  m ,  ) , (A-8.4) X

X

n

with the decomposition coefficients * uˆ( xn ,  m ,  )  u( x,  ), bmn ( x )   u( x,  )  bmn ( x ) dx , where

,

(A-8.5)

stands for inner product, “*” stands for complex conjugate, and bmn ( x ) are the dual

GD frame atoms,

bmn ( x )  g ( x  xn )eim x , (A-8.6) with g ( x ) as the dual-window function of g ( x ) . The dual-window function can be calculated by pseudo-inversion of the original window function (Qian and Chen, 1996; Mallat, 1999; Wu and Chen, 2001). From equation (A-8.4), we can see that the local plane wave is a weighted average of the windowed beamlets with same wavenumber from neighboring windows. The total space domain wavefield can be written as u ( x,  )   u ( x, m ,  ) . (A-8.7) m

Appendix B: Local Exponential Frame (Lef) Beamlet Decomposition The LEF (e.g., Daubechies et al., 1991; Auscher, 1994; Mao and Wu 2007) beamlet decomposition can eliminate the directional symmetry in the local cosine basis (LCB) beamlets. Same as the GDF beamlet, it can decompose the wavefield into LWD with uniquely defined direction information. The LEF beamlets for the beamlet decomposition in equation (8.1) are bell function modulated exponential harmonics,

97

bmn ( x ) 

2 Bn ( x ) exp  i m ( x  xn )  , Ln

(B-8.1)

where Bn ( x ) is a smooth bell (window) function (fore detail, see Wu et al., 2008), Ln  xn 1  xn is the nominal length of the window, and  m    m  12  Ln

( m   M  1,  ,  1, 0, 1,  , M ). Above LEF beamlets can be separated into the right- and () ( ) ( x ) and bmn ( x) , left-propagating beamlets bmn  () bmn ( x )    b(  ) ( x )   mn 

2 (c) (s) Bn ( x ) exp  i m ( x  xn )   bmn ( x )  ibmn ( x) Ln 2 (c) (s) Bn ( x ) exp  im ( x  xn )   bmn ( x )  ibmn ( x) Ln

,

(B-8.2)

 (c) 2 Bn ( x ) cos m ( x  xn )  bmn ( x )  Ln  , (B-8.2)  2 b( s ) ( x )  Bn ( x ) sin  m ( x  xn )   mn Ln  (c) (s) where m  1, 2,  , M now, and bmn ( x ), bmn ( x ) are the LCB and LSB (local sine basis) beamlets respectively. Substituting equation (B-8.2) into (8.2), we can obtain the partially reconstructed right- and left-propagating local plane waves u (  ) ( x,  m ,  ) and u (  ) ( x,  m ,  ) ,  () 2 Bn ( x )uˆ (  ) ( xn ,  m ,  ) exp( im xn ) u ( x, m ,  )  exp( im x ) Ln n  ,  2 u (  ) ( x,  ,  )  exp( i x ) () n L Bn ( x)uˆ ( xn , m ,  ) exp( im xn ) m m  n  with the decomposition coefficients () uˆ (  ) ( xn , m ,  )  u( x,  ), bmn ( x)  .  () () uˆ ( xn ,  m ,  )  u( x,  ), bmn ( x ) The space-domain total wavefield can be written as u ( x,  )   u (  ) ( x,  m ,  )   u (  ) ( x,  m ,  ) . m

(B-8.4)

(B-8.5)

(B-8.6)

m

Using the relation between the LEF and the LCB/LSB beamlets in (B-8.2), the LEF () () decomposition coefficients uˆmn and uˆmn in equation (B-5) can be represented with the LCB and X

(c) mn

LSB decomposition coefficients uˆ

X

(s) mn

and uˆ , (c) (s)  (  ) uˆmn  iuˆmn ˆ umn  4 .  (c) (s) ˆ ˆ  u iu mn uˆ (  )  mn  mn 4

98

(B-8.7)

The LCB and LSB decomposition coefficients can be obtained by fast algorithms, such as the one introduced by Malvar (1992) (see also Mallat, 1999), which makes the LEF decomposition coefficients can be obtained very efficiently.

99

9. A DEPTH MIGRATION METHOD BASED ON THE FULLWAVE REVERSE-TIME CALCULATION AND LOCAL ONEWAY PROPAGATION 9.1. Background The full-wave reverse-time migration method does not have the angle limitation of one-way propagators and can be used to image structures in complex areas especially for steeply dipping structures where turning waves are involved. Although the reverse-time migration is a relatively expensive method, the development of computer hardware and software (remarkably the popularities of PC clusters and Linux operating system) plus the ever increasing interests to explore complex structures in challenging regions stimulated recent interest to this method (e.g., Wu, et al., 1996; Mufti, et al., 1996; Sun and McMechanz, 2001; Biondi and Shan, 2002; Yoon, et al., 2003, 2004; Mulder and Plessix, 2004; Fletcher, et al., 2005; Wang, et al., 2005; Hou and Marfurt, 2005). One of the major problems that affect the application of reverse-time migration is that the wideangle capability of full-wave propagator can generate spurious cross correlations from diving waves and back-scattered waves, and causes serious artifacts (Yoon, et al., 2004; Mulder and Plessix, 2004; Fletcher, et al., 2005; Wang, et al., 2005). Special techniques such as velocity model smoothing or interface impedance matching have been introduced to solve this problem. These approaches can improve the situation in poststack migration where only one wavefield is involved, but the result is unsatisfied when applied to the prestack migration where cross correlation happens between the source and receiver waves. Recently, several approaches have been proposed to solve this problem. Mulder and Plessix (2004) suggested that these artifacts can be removed by blanking the data, high-pass filtering or iterative migration. Yoon, et al. (2004) proposed to use an angle related weighting function calculated from Poynting vectors of source and receiver waves to pick the right events for imaging. Fletcher, et al. (2005) tried to eliminate the artifacts by introducing a directional damping term in areas where unwanted reflections occur. Wang et al. (2005) proposed a pseudospace method by using a space distortion to transform a heterogeneous velocity model into a homogeneous model where no reflections were generated. In this paper, we propose to use the local one-way propagator as the filter to separate the fullwave equation generated wavefield into directional waves. Then use directional wavefields for imaging. By making proper combinations of these directional waves, the artifacts caused by undesired waves can be avoided. Another advantage of this method is that it only needs to store part of the wavefield. The size of the file containing the wavefield can be reduced by about an order-of-magnitude which dramatically reduces the effort in input/output, network traffic and data management.

9.2. Methodology

100

The wavefield in a closed region without internal source can be calculated from its boundary values using the Rayleigh integral G  r, ξ  u  ξ     G  r, ξ  u  r   S  u  ξ  ds , n n  

(9.1)

where u  r  is the wavefield, G  r, ξ  is the Green function, S is the boundary, r is the position of an internal point, ξ is the position on the boundary, and  n is the derivative along the boundary normal. From the wavefield calculated using the full-wave propagator, choose the wavefield along certain boundaries, the wavefield in the entire region can be recovered from these boundary values. If we choose the boundaries along different directions and use one-way propagator as the Green function in equation (9.1), we can obtain waves propagating along different directions. For convenience, we use the wavefield on horizontal lines (or planes for 3D case) for wavefield extrapolation in near vertical direction and use wavefield on vertical lines for near horizontal extrapolation. By choosing different normal directions of the boundary and using forward or backword one-way propagators, we can extrapolate wavefield in forward or backword directions.

Figure 9.1. A simple two-layer velocity model. The horizontal and vertical lines indicate the locations where wavefields are needed for one-way extrapolations. Taking the horizontal line as an example, the wavefield u  ξ  and u  ξ  z are obtained from the full-wave calculation. Substituting these values into equation (9.1) and using one-way Green function G  z  r, ξ  and G  z  r, ξ  z , where superscript  z denotes the one-way propagator along positive z direction, we can calculate wavefield below that horizontal line. Similarly, using Green function G  z  r, ξ  and G  z  r, ξ  z , we can calculate wavefield above the horizontal line and along  z direction. Using wavefield along horizontal lines at different depth levels, we can reconstruct wavefields along  z and  z directions between these lines. A similar method can be used to process the near horizontally propagating waves. In these cases, the one-way propagator only deals with narrow-angle propagation and the wavefield is extrapolated only for a very short distance. The accuracy of the one-way propagator will not seriously affect the result.

101

Figure 9.2. Wavefield snapshots calculated using full-wave finite-difference method. (a) and (b) are source wavefields at 0.8 s and 1.2 s. (c) and (d) are receiver wavefields at T  0.8s and T  1.2s .

9.3. Numerical Examples To show how this system works, we first conduct the migration for a simple two-layer model. Figure 9.1 shows the velocity model. The size of the model is 451150 with dx  dz  24 m , and the velocities are 3.0 km/s and 4.0 km/s for upper and lower layers, respectively. The horizontal and vertical lines in the velocity model indicate the location where the wavefield is picked. Both the synthetic data and the reverse-time migration are calculated using a fourth-order full-wave finitedifference code. To focus our attention to the effect of back-scattered waves, We mute the first arrival from the shot record. Shown in Figure 9.2 is migrated wavefield using the full-wave method, with 2a and 2b are source waves at 0.8 s and 1.2 s, and 2c and 2d are receiver waves at T  0.8s and T  1.2s , where T is the maximum recording time. We then output the wavefield on horizontal lines at depths for every 15 dz and use a one-way propagator (Xie and Wu, 1998, 2005) to extrapolate the wavefield from each depth to the next depth (15 steps). Figures 9.3a and 3b are the reconstructed wavefield us z  r  , and 3c and 3d are the reconstructed wavefield u g z  r  , where

102

subscripts s and g denote source and receiver waves. Compared with Figure 9.2, the undesired reflections have been eliminated. The image can be calculated from directional wavefields. To avoid artifacts from cross correlations between diving waves and back-scattered waves, we choose wave pairs related to the true reflections, for example us z  r  u g z  r  , us z  r  u g z  r  , us x  r  u g x  r  or us x  r  u g x  r  , to form the image. The final image can be obtained by summing up all partial images. Shown in Figure 9.4 are reverse-time images for the two-layer model, where 4a is calculated using the conventional zero-lag cross correlation image condition T (9.2) I  r   0 u s  r, t u g  r , T  t  dt , where I  r  is the image,

us and u g

are source and receiver waves. Figure 9.4b is calculated using

the image condition given by Yoon et al. (2004). In Figure 9.4c, the image is calculated from directional waves us z  r, t  ug z  r, T  t  reconstructed using local one-way propagator. In Figure 9.4a, there are artifacts shown above the interface, while in 9.4b and 9.4c, they have been properly removed. As second example, we calculate the image of a simple salt dome model. Figure 9.5 shows the velocity model overlapped with the wavefield snapshot from the finite-difference forward modeling. The velocity model composed of a background with a vertical velocity gradient and a salt dome with steep flanks. The size of the model is the same as the two-layer model and major phases from different part of the structure are indicated with rays. Figure 9.6 shows the wavefield from finite-difference forward modeling. In the successive time slices, we can see the large angle reflections developed from the salt flank. Even for such a simple velocity model, the reflections are rather complex. To properly image the steep flank of the salt dome, wide-angle turning waves must be properly handled. Figure 9.7 shows the images from depth migrations overlapped with the velocity model. A single source is used to illuminate the model. Figure 9.7a is calculated using the conventional zero-lag cross correlation imaging condition. Figure 9.7b is calculated using the imaging condition of Yoon et al. (2004). Figure 9.7c is calculated from the waves us z  r  and u g z  r  reconstructed using the local one-way propagator, and 9.7d is calculated using us x  r  and u g x  r  . Note that horizontal structure is better imaged in 9.7c and vertical structure is better imaged in 9.7d. The artifacts shown in 9.7a have been removed in 9.7b to 9.7d.

103

Figure 9.3. Wavefield snapshots reconstrcted using the one-way propagator. (a) and (b) are wavefields along  z and at 0.8 s and 1.2 s. (c) and (d) are wavefields along  z and at T  0.8s and T  1.2s . Note that compared with the finite-difference result shown in Figure 9.2, the undesired reflections have been eliminated.

104

Figure 9.4. Reverse-time images for the two-layer model, (a) using conventional zero-lag image condition, (b) using the image condition given by Yoon et al. (2004), and (c) using full-wave propagator plus local one-way propagator.

Figure 9.5. A simple salt dome model with steep flanks. Overlapped is the snapshot generated by finite-difference forward modeling. The major reflection events are marked with rays.

9.4. Discussions and Conclusion 105

We proposed a depth migration method based on the full-wave reverse-time calculation and local one-way propagation. Both the source and receiver waves are extrapolated using the full-wave finite-difference propagator. There is no angle limitation for these wavefields. Then the local one-way propagator is used to separate the waves into different directions. Finally, the image is calculated using directional wavefields. Using the one-way propagator at the imaging stage has the following advantages. (i) It eliminates undesired waves and removes the artifacts. (ii) The input/output file size can be substantially reduced (about an order-of-magnitude smaller than if the entire wavefield is dumped out). The additional computation for the one-way propagation can be compensated by the time saved from the I/O and network traffic. The reduced file size makes the management of the data easier and provides the possibility for repeated later studies.

Figure 9.6. Snapshot of finite-difference forward modeling. Note the complex large angle reflections from the salt flank. The downward arrival at the later time is the artifacts from the upper boundary. The result presented in this paper is preliminary. Further investigation is needed to better reconstruct the wavefield in complex regions using one-way propagator. More sophisticated method such as the angle domain image analysis (Wu, et al. 2004, and Xie, et al., 2005) may be required to combine images from multi-directional wavefileds.

106

Figure 9.7. Depth migration image for a simple sand dome model, (a) using the conventional zero-lag imaging condition, (b) using the image condition of Yoon et al. (2004), (c) using us z  r  u g z  r  u x r ux r and (d) using s   g   .

107

PART III. IMAGING AND INVERSION

108

10. GENERALIZED DIFFRACTION TOMOGRAPHY IN HETEROGENEOUS MEDIA WITH FINITE DATA APERTURE 10.1. Backpropagation and the Imaging Principle The traditional diffraction-tomography with a homogeneous reference model applies a filtered backpropagation to the scattered field (data) for the reconstruction of the parameter perturbations in the model space (Devaney, 1982, 1984; Beylkin et al., 1986; Wu and Toksöz, 1987, Harris, 1987; Pratt and Worthington, 1988: Wu et al., 1994). The backpropagation process is a doublefocusing operation, which focus both the waves from the source array and the receiver array to the image point. Therefore it is equivalent to applying the imaging condition at the image point in the migration process (Claerbout, 1985). For homogeneous background with infinite acquisition aperture, the propagation of different plane-wave components is independent of each other, and the deconvolution filter (decon-filter) remains unchanged across the whole image space. Therefore the order of propagation and filtering in wavenumber domain is interchangeable and the filtering is applied to the scattering data before the backpropagation. For heterogeneous backgrounds and/or finite data aperture, this symmetry is broken and the filter is location dependent. For inhomogeneous background, different plane wave components in the data will be distorted and cross-coupled with each other. The simple relation between the data spectra and the medium perturbation spectra no longer exists. In order to overcome the limitation of the classic diffraction tomography, some spatial domain formulation (Woodward, 1989) or waveform inversion approach (Pratt and Worthington, 1990; Zhou et al., 1995, 1997; Xu et al., 1995; Pratt, 1999; Min and Shin, 2006; Sheng et al., 2006) were developed for the case of inhomogeneous background media. However, in these approaches, the intuitive spectral inversion and its efficiency had been lost. In this paper, I will generalize the diffraction tomography to the case of inhomogeneous background with one-way propagators (Green’s function) and resolving kernels in the local angle domain. The process is to first perform the backpropagation of the scattered field for removing the background effect and then apply the space-dependent decon-filter in the local wavenumber domain. We will see that when the aperture is infinite in a homogeneous background, the process reduces to the filtered backpropagation of the classic diffraction tomography. The parameters to be inverted are volume perturbations of unknowns with respect to a given reference medium. The Green’s function is the impulse response in the reference medium. In the Born model, the scatterings of each volume elements are independent from each other and no interaction between the elements is taken into consideration. For the case of scalar media (acoustic media with constant density), the scattered (pressure) field under the Born approximation is psc ( xg , xs ) = -k 2 ò GM ( x, xs )GM ( xg , x )O ( x )dx (10.1) V

where O (x) = 1 -

C 02 C 2 (x)

109

(10.2)

is the velocity perturbation function of the medium, which is called the object function in classic diffraction tomography (Devaney, 1982, 1984; Wu and Toksöz, 1987), and GM is the scalar wave Green’s function for the modeling process in the inhomogeneous background media. The exact GM should be a full wave Green’s function. However, it can be replaced by different approximations in different inversion schemes. The Born approximation itself is only valid for weak scattering, i.e. when the parameter perturbations are small and the integration volume (heterogeneity extension) is small. However, since GM can be updated and set to be close to the Green’s function for the real media, (10.1) can be considered as a local Born approximation which is valid for large volume scattering problems.

10.2. Imaging Condition in the Local Angle-Domain and the Local Image Matrix The standard imaging condition in the space domain is in the form of cross-correlation, L(w, x) = 2 ò

A(xs )

dxsGI* (w, x,; xs ) ⋅ ò

A(x g )

dx g

¶GI* (w, x,; x g ) sc p (w, x g ; xs ) ¶z

(10.3)

where GI is Green’s function used in the imaging process, which could be different from the Green’s function of forward modeling; “ * ” stands for complex conjugate. The inner integral is a back propagation Rayleigh integral and A(x g ) is the spatial receiver aperture for the given source. The outer integral is the summation over all the sources and A(xs ) is source aperture. This is the double focusing for a single-frequency wavefield. In migration imaging, the final image is obtained by summing up images of all the frequencies, a double-focusing in time-domain. In diffraction tomography, both single-frequency tomography and multi-frequency tomography can be implemented. In order to be symmetric for the source array focusing and receiver array focusing and compare later with the classic diffraction tomography, we slightly modify the imaging condition (10.3) into L(w, x) = 4 ò dxs As

¶GI* (w, x,; x g ) sc ¶GI* (w, x,; xs ) ⋅ ò dx g p (w, x g ; xs ) Ag ¶z ¶z

(10.4)

In this imaging condition, we propagate the source field as if from boundary elements instead of volume elements (isotropic source). Later we will see the symmetry and convenience of this modification. In order to obtain the local angle-spectra of an image field, imaging condition has been extended from space domain to the space-angle-domain (beamlet domain) (Wu and Chen, 2002, 2004). Then the image function is no longer a scalar value and becomes a matrix: LIM (local image matrix) L(x, ks , kg ) , where x = (x, z ) is the position vector at depth z; ks and kg are the source and receiving wavenumbers, respectively. To simplify the notations, we omit the frequency dependence in the LIM and the following derivations. When necessary, we will specify the frequency dependence again explicitly. Note that the source direction is defined as the direction from the image point to the source on the surface, and is opposite to the incident direction. The new imaging condition (for a single frequency) in the space-angle domain is 110

L(x, ks , kg ) = 2ò

A(xs )

dxsGI* (x, ks ; xs ) ⋅ ò

A(x g )

dx g

¶GI* (x, kg ; x g ) sc p (x g ; xs ) ¶z

(10.5)

where GI (x, ks ; xs ) is the incident beamlet at the image point generated by a point source at xs on the surface, and GI (x, kg ; x g ) , the outgoing (scattered) beamlet at the imaging point received by a point receiver at x g on the surface. If the Green’s functions used in the inversion are exactly the same Green’s functions as those in the forward modeling (in the acquisition process), or their phase information are exact, the LIM will be a real number matrix. However, due to the approximations in Green’s function and velocity model inaccuracy, the LIM in general is a complex matrix: L(x, ks , kg ) = [ AL exp(iFL ) ] , (10.6) where AL = AL (x, ks , kg ) is the amplitude of matrix element and FL = F(x, ks , kg ) is the phase of matrix element. For an inaccurate background velocity model, the phase FL is different for different pairs of ( ks , kg ), resulting in incomplete focusing. These phase residues can cause the blurriness of the final image and produce artifacts and coherent noises. However, under the Born approximation, these effects are assumed small and negligible. Furthermore, the image phaseresidue matrix can also be used to update the velocity model during an iterative process.

10.3. Deconvolution Filtering in Local Wavenumber Domain and the Generalized Diffraction Tomography The modeling and imaging processes can be written into operator forms: (10.7) I(x) = B(x | w, xs , x g )U(w, xs , x g ) (10.8) where F is the acquisition (modeling) operator which maps the model S into the data set U; while B is the imaging operator which invert the data U into the image I, x0 is the scattering point in the model space, and x is the location of the image point in the imaging process. Substitute (10.7) into (10.8) we have (10.9) I(x) = B(x | w, xs , x g )F(w, x g , xs | x0 )S(x0 ) . The resolution operator (matrix) is obtained as R(x, x 0 ) = B(x | w, xs , x g )F(w, x g , xs | x 0 ) (10.10) If the imaging operator is exactly the inverse operator of the modeling operator, the resolution operator will be an identity operator. For most cases, the resolution matrix is not an identity matrix and the spreading of the matrix elements along the diagonal give some quantitative measure of the parameter resolution of the imaging/inversion. If the resolution matrix can be calculated, the true model can be obtained by deconvoluting the image field with the resolution matrix (with some regularization). However, the calculation of resolution operator and volume deconvolution in space-domain is intractable. We will formulate the process of decon-filtering in the local wavenumber domain, i.e. the beamlet domain. U(w, x g , xs ) = F(w, x g , xs | x0 )S(x0 )

In the following we will derive the resolution operator for the imaging process based on the Born model for volume perturbations. After decomposing the Green’s functions into the local wavenumber-domain at a level z near the image point x ' = (x ', z ') , and substituting into (10.4), the image matrix at x ' becomes

111

L(x ', ks , kg ) = - exp{-i(kg + ks ) ⋅ r '}4Vs Vg ò dxsGI* (x, ks ; xs ) As

⋅ò

Ag

dx gGI* (x, kg ; x g )psc (x g ; xs )

,

(10.11)

where r ' = (x ', z '- z ) and ks = (xs , Vs ) , kg = (xg , Vg ) , with x and V as the corresponding local horizontal and vertical wavenumbers respectively. From equation (10.11), we see that from a pair of local incident-scattering angles we can only detect the local spectral component at K = ks + kg = k ^ . Set z ' = z , equation (10.11) can be written into another form L(x ', K = k ^ ) = -4Vs V g exp{-i(xg + xs ) ⋅ x '}ò dxsGI* (x ', ks ; xs ) As

⋅ò dx gGI* (x ', kg ; x g )psc (x g ; xs )

.

(10.12)

Ag

Figure 10.1 shows The definition of local wavenumbers and the domain of spectral coverage (domain of integration) for the Born model. With multi-frequency imaging, the spectral coverage can be expanded and the local image matrix becomes I (x ', K) = ò d wL(w, x ', K = k ^ ) . (10.13) Af

Figure 10.1: The definition of local wavenumbers and the domain of spectral coverage for the Born model The angular spectra of local images are usually highly non-uniform due to the acquisition geometry, limited data aperture and the influence of overburden structures. In order to reconstruct the true parameter of the perturbation function (the object function), we need to deconvolute the local image field with the resolving kernel. The kernels for the resolution operator (resolving kernel) (Backus and Gilbert, 1970; Tarantola, 1987, 2005) is obtained based on (10.10), (10.4) and (10.1) by substituting the forward modeling (the Born model) (10.1) into (10.10):

112

Â(x; x 0 ) = -4 ò d wk 2 ò dxs As

¶GI * (w, x; xs ) GM (w, x 0 ; xs ) ¶z

*

¶GI (w, x; x g ) GM (w, x g ; x0 ) òAg dx g ¶z

(10.14)

For detailed derivation of the resolution operator and its kernel, see Wu et al. (2006). The resolving kernel is also called the point spreading function (PSF) for the imaging system which includes both the acquisition (modeling) and inversion (imaging) processes. We will formulate the process in the local wavenumber domain, i.e. the beamlet domain. We perform local 3D Fourier transform (local wavenumber-domain decomposition) on Â(x, x0 ) with coordinate center at x 0 , resulting in the expression of PSF in local wavenumber domain (Wu et al., 2006) Â(K = k ^ ; x 0 ) = -2ò d wk 2 ò dxs Af

As

¶GI * (w, x0 , ks ; xs ) GM (w, x 0 ; xs )BI (x ', kg ) ¶z

(10.15)

where BI (x ', kg ) is the backpropagation integral of the modeling Green’s function, BI (x ', kg ) = -2 ò dx g Ag

¶GI* (w, x ', kg ; x g ) GM (w, x '; x g ) . ¶z

(10.16)

Beamlet decomposition (localized Fourier decomposition) decomposes the field of a Green’s function, a spherical wave front, into beamlets. Each beamlet has a center location x 0 with certain width, and a lobe in the specified direction. Because the existence of the lobe, the directional resolution is not as sharp as in the global Fourier transform. The spatial localization of the beamlet is not a point, but with certain width. The compromise between the spatial and directional localizations is governed by the uncertainty principle. ¶GI * (w, x 0 , ks ; xs ) GM (w, x 0 ; xs ) can be considered as a point source at ¶z to xs , and then backpropagated (phase conjugate) again to x 0 , with the local

In (10.15) and (10.16),

x 0 radiating wavefield plane wave ks . If we use the high-frequency asymptotic approximation of Green’s function, then the incident wave direction will be determined by the ray direction. For the more general case of wave propagation, the local plane wave decomposition is more appropriate (Wu and Chen, 2002;

¶GI * (w, x ', kg ; x g ) GM (w, x '; x g ) corresponds to ¶z a point source at x ' , radiating wavefield to x g , and then backpropagated to x ' , with the local plane

Xie and Wu, 2002; Wu et al. 2004). In a similar way

wave kg . The acquisition aperture effects are included in the integrals with apertures As ,

Ag

and Af .

The image field I is seen from (10.9) a convolution of the resolution operator R with the model S. Therefore, the true object function is obtained by deconvoluting the image with the resolving kernel. In the local wavenumber domain, the decon-filter is a division of (10.13) by (10.15): O(x, K) = I (x, k ^ )/ Â(k^ ; x) (10.17) The local image in the space-domain can be reconstructed by the inverse local Fourier transform, O(x ') = òò dKe iKx x 'O(x ', K) . (10.18) The above approach is to perform the decon-filtering to the final multi-frequency image. An alternative approach is to filter the images in the local wavenumber domain for each single

113

frequency, and obtain the final object spectrum by averaging over the single-frequency object spectra. For a single-frequency, the reconstruction (10.17) becomes

òò dKe O(x ', k^ ) = òò dKx dKze O(x ', k^ ) , = òò d xsd xge i(x + x )x 'J (K x , K z ; xs , xg )O(x ', k ^ ) = òò d xsd xge i(x + x )x '  (xs , xg )L(x ', xs , xg )

O(x ') =

iKx x '

iKx x '

s

g

s

g

(10.19)

which is seen to be a backpropagation plus filtering process. In (10.19) J (Kx , Kz ; xs , xg ) is the Jacobien of the coordinate transform from ( xs , xg ) to ( K x , K z ). The filter  (xs , xg ) =

J (K x , K z ; xs , xg )

Â(k ^ ; x ')

=

xs Vg - xg Vs Vs V g Â(xs , xg ; x ')

(10.20)

is the filter (tomographic filter) which can be applied directly to the image matrix in the local wavenumber domain to get the object function. In the following, we will further simplify the filter function by considering the behavior of resolving kernel in the beamlet domain. In (10.15) and (10.16), we can also put the beamlet decomposition to the modeling Green’s functions GM (w, x 0 ; xs ) and GM (w, x '; x g ) instead of the imaging Green’s functions and consider the integrals as limited-aperture backpropagation integrals. Then we can define the modified resolving kernel for a single frequency as Vs Vg Â(xs , xg ; x ') = 4k 2 ò dxs As

¶GI * (w, x0 ; xs ) ¶GM (w, x0 , ks ; xs ) ¶z ¶z

¶GI* (w, x '; x g ) ¶GM (w, x ', kg ; x g ) dx g òAg ¶z ¶z

(10.21)

and the tomographic filtering becomes a more stable process without singularities in spectral domain:  (xs , xg ) =

xs Vg - xg Vs k2

/ A(xs , xg )

(10.22)

where A(xs , xg ) is factor quantifying the acquisition aperture effect in the given inhomogeneous background: A(xs , xg ) = 4 ò dxs As

¶GI * (w, x0 ; xs ) ¶GM (w, x0 , ks ; xs ) ¶z ¶z

¶GI* (w, x '; x g ) ¶GM (w, x ', kg ; x g ) òAg dx g ¶z ¶z

(10.23)

In case we know the velocity model of the medium to be imaged, we can simply the calculation of (10.23) by using the efficacy matrix in the illumination analysis (Wu and Chen, 2002, 2006). The back propagation integral in (10.23) represents a refocusing process of the forward propagated beamlet (spreading) wavefront with limited wave front aperture. For large enough receiver (or source) aperture, we consider that all the energy collected by the receiver array can be refocused to the lobe of the beamlet. Then energy conservation exists GI* ( , x '; xg ) GM ( , x ', k g ; xg ) 2 G ( , x ', k g ; xg ) 2 (10.24) |  dxg |   dxg | M | Ag Ag z z z

114

Apply the same approximation for the source array backpropagation integral, the aperture factor (10.23) can be calculated by a similar procedure for the AAE (acquisition aperture efficacy): A(xs , xg )

2

=

òA dxs s

|

¶GM (w, x ', xs ; xs ) 2 | ¶z

òA dx g

¶GM (w, x ', xg ; x g ) 2 | ¶z

|

g

(10.25)

10.4. Relation with the Classic Diffraction Tomography in A Homogeneous Background In a homogeneous media with infinite aperture of acquisition, the local wavenumber becomes global and the beamlet decomposition is replaced by the global Fourier transform. Then the Green’s function in (10.12) becomes G I ( ks ; x s ) =

i exp ( iVs z ) exp ( i xs xs ) . Vs 2

(10.26)

Substituting the Green’s function into (10.12), resulting in L(z ', K = k ^ ) = e -i(Vs + Vg )z ' ò dxse -ixs xs ò dx ge -i xg x g psc (x g ; xs ) = e -i(Vs + Vg )z ' psc (xg ; xs )

.

(10.27)

This is a simple backpropagation of the scattered field in the wavenumber domain. In a similar way, the decon-filter in (10.15) for a single-frequency becomes Â(k ^ ; x 0 ) = k 2 ò dxse -i xs xs GM (w, x 0 ; xs )ò dx ge -i xg xg GM (w, x g ; x 0 ) =

k2 Vs V g

.

(10.28)

Figure 10.2 shows the inverse of the resolving kernel. The tomographic inversion in wavenumber domain becomes O(k ^ ) = L(k ^ )/ Â(k ^ ) = e -i(Vs + Vg )z '

Vs V g k2

psc (xg ; xs )

(10.29)

Transforming back to the space-domain, we get the tomographic image of the object function: O(x ') = =

òò dKO(k^ ) = òò dKx dKzO(k^ )

òò d xsd xgJ (Kx , Kz ; xs , xg )e

-i (Vs + V g )z ' Vs V g 2

=

òò d xsd xge

i (xs + xg )x ' -i (Vs + V g )z '

=

òò d xsd xge

i (xs + xg )x ' -i (Vs + V g )z '

e e

psc (xg ; xs ) k . xs V g - xg Vs sc p (xg ; xs ) k2

(10.30)

 (xs , xg )psc (xg ; xs )

where  (xs , xg ) is the tomographic filter function  (xs , xg ) =

J (K x , K z ; xs , xg ) Â(k ^ ; x ')

=

xs V g - xg Vs k2

(10.31)

as shown in Figure 10.3. Note that with infinite aperture in homogeneous media, the general equation of tomographic filter in (10.22) converges to (10.31). Since for normalized Green’s function in beamlet domain, the integral in (10.25) approaches unity. 115

The above formulation is exactly the filtered backpropagation of classic diffraction tomography in homogeneous backgrounds for surface reflection data with an infinite aperture (Devaney, 1982, 1984; Wu and Toksöz, 1987). In the more general case of heterogeneous backgrounds, as we proposed above, the process is a backpropagation plus filtering (deconvolution) in the local angle-domain.

Figure 10.2: The spectrum weighting function, which is the inverse of resolving kernel.

Figure 10.3. The tomographic filter function in the local wavenumber domain

10.5. Conclusion Tomographic inversions in generally inhomogeneous media can be summarized as a process of backpropagation plus filtering in the local angle domain. The backpropagation is a doubly focusing process, similar to the imaging principle in migration/imaging with extended imaging conditions in the local wavenumber domain. The filtering is a deconvolution in the local angle domain. The generalization of the classic diffraction tomography enables the application of the method to more general cases with arbitrary acquisition geometry. Since the inhomogeneous background can be set to be close to the real media, the applicability of the Born or Rytov approximations are also extended.

116

REFERENCES Aki, K. and Richards, P.G., 1980, Quantitative Seismology: W.H. Freeman and Company, San Francisco. Aminzadeh, F., N. Burkhard, L. Nicoletis, F. Rocca, and K. Wyatt, 1994, SEG/EAGE 3-D modeling project: 2nd update: The Leading Edge, 13, 949–952. Aminzadeh, F., N. Burkhard, T. Kunz, L. Nicoletis, and F. Rocca, 1995, 3-D modeling project: 3rd report: The Leading Edge, 14, 125–128. Aprea, C. S. Hildebrand, M. Fehler, L. Steck, W. Baldridge, P. Roberts, C. Thurber, and W. Lutter, 2002, Three-dimensional Kirchhoff migration: Imaging of the Jemez volcanic field using teleseismic data, J. Geophys. Res. doi:10.1029/2000JB000097. Audebert, F., P. Froidevaux, H. Rakotoarisod, and J. Svay-Lucas, 2002, Insights into migration in the angle domain: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1188-1192. Auscher, P., 1994, Remarks on the local Fourier bases, in J.J. Benedetto and M.W. Frazier, eds., Wavelets, Mathematics and applications: CRC Press, 203-218. Backus, G. and F. Gilbert, 1968, The resolving power of gross Earth data, Geophys. J. Royal Astr. Soc. 13, 247-276. Bear, G., C. Liu, R. Lu, D.Willen and I. Watson, 2000, The construction of subsurface illumination and amplitude maps via ray tracing: The Leading Edge, 19, 726–728. Berkhout, A.J., 1997, Pushing the limits of seismic imaging, Part I: Prestack migration in terms of double dynamic focusing: Geophysics, 62, 937-954. Berkhout, A. J., Ongkiehongz, L., Volker, A. W. F., and Blacquière, G., 2001, Comprehensive assessment of seismic acquisition geometries by focal beams—Part I: Theoretical considerations, Geophysics, 66, 911-917. Beylkin G., Oristaglio M. and Miller D. 1985, Spatial resolution of migration algorithms, in Proceedings of the 14th International Symposium on Acoustic Imaging, pp. 155-167. Plenum Press. Beylkin, G., 1985, Imaging of discontinuities in the inverse scattering problem by inversion of a causal generalized Radon transform, J. Math. Phys., 26, 99-108. Beylkin, G., Otistaglio, M. and Miller, D., 1986, Spatial resolution of migration algorithms, Acoustic imaging, v. 14, 155-168. Billette, F. J., and S. Brandsberg-Dahl, 2005, The 2004 BP velocity benchmark: 67th Meeting, EAGE, Expanded Abstracts, B035. Biondi, B. and G. Shan, 2002, Prestack imaging of overturned reflections by reverse time migration: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1284-1287. Biondi, B., 2007, Angle-domain common-image gathers from anisotropic migration: Geophysics, 72, S81-S91.

117

Biondi, B., and W. Symes, 2004, Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging: Geophysics, 69, 1283–1298. Bleistein, N., J. K. Cohen, and F. G. Hagin, 1987, Two and one-half dimensional Born inversion with an arbitrary reference: Geophysics, 52, 26–36. Bleistein, M., J.K. Cohen and J.W. Stockwell, Jr, 2001, Mathematics of multidimensional seismic imaging, migration and inversion, Springer. Blomgren, P, G. Papanicolaou and H Zhao, 2002 Super-resolution in time-reversal acoustics, J. Acoust. Soc. Am. 111, 230-248. Brandsberg-Dahl, S., M.V. de Hoop, and B. Ursin, 2003, Focusing in dip and AVA compensation on scattering angle/azimuth common image gathers: Geophysics, 68, 232254. Cao, J., and R.S. Wu, 2005, Influence of Propagator and Acquisition Aperture on Image Amplitude: 75th Annual International Meeting, SEG, Expanded Abstracts, 1946-1949. Cao, J., and R.S. Wu, 2006, Study of the influence of propagator amplitude correction on image amplitude using beamlet propagator with local WKBJ approximation: 76th Annual International Meeting, SEG, Expanded Abstracts, 2499-2503. Cao, J., and R.S. Wu, 2008, Amplitude compensation for one-way wave propagators in inhomogeneous media and its application to seismic imaging: Communications in computational physics, 3, 203-221. Cao, J., and R.S. Wu, 2008, Fast acquisition aperture correction by beamlet decomposition: 78th Annual International Meeting, SEG, Expanded Abstracts, submitted. Cao, J., and R.S. Wu, 2008, Local-angle domain illumination for full-wave propagators, 78th Annual International Meeting, SEG, Expanded Abstracts, submitted. Chen, J. and Gerard T. Schuster, G., 1999, Resolution limits of migrated images, Geophysics 64, 1046-1053. Chen, J. and Gerard T. Schuster, G., 2001, On "Resolution limits of migrated images" (J. Chen and G.T. Schuster, Geophysics, 64, 1046-1053). Reply Geophys. 66, 693-694 . Chen, L., R.S. Wu, and Y. Chen, 2006, Target-oriented beamlet migration based on GaborDaubechies frame decomposition: Geophysics, 71, S37–S52. Claerbout, J. F., 1971, Toward a unified theory of reflector mapping: Geophysics, 36, 467–481. Claerbout, J.F., 1985, Imaging the earth interior: Blackwell Scientific Publications. Clayton, R.W., and R.H. Stolt, 1981, A Born-WKBJ inversion method for acoustic reflection data: Geophysics, 46, 1559-1567. Coifman, R.R. and Y. Meyer, 1991, Remarques sur l’analyse de Fourier a fenetre, Comptes Rendus de l’Academie des Sciences: Paris, Serie I 312, 259-261. Collins, M.D., 1989, Applications and time-domain solution of higher-order parabolic equations in underwater acoustics, Journal of the Acoustical Society of America, 86, 1097-1102.

118

Collins, M. D., and E. K. Westwood, 1991, A higher-order energy-conserving parabolic equation for range-dependent ocean depth, sound speed, and density: Journal of the Acoustical Society of America, 89, 1068-1175. Daubechies, I., S. Jaffard, and J.-L. Journé, 1991, A simple Wilson orthonormal basis with exponential decay: SIAM Journal of Mathematical Analysis, 22, 554-573. de Bruin, C.G.M., C.P.A. Wapenaar, and A.J. Berkhout, 1990, Angle-dependent reflectivity by means of prestack migration: Geophysics, 55, 1223-1234. De Hoop, M., J. Le Rousseau, and R. S. Wu, 2000, Generalization of the phase-screen approximation for the scattering of acoustic waves: Wave Motion, 31, 43-70. Devaney, A. J., 1981, Inverse-scattering theory within the Rytov approximation: Optics Letters, 6, 374-376. Devaney, A. J., 1982, A filtered back propagation algorithm for diffraction tomography: Ultrasonic Imag., 4, 336-350. Devaney, A. J., 1984, Geophysical diffraction tomography: Trans., Inst. Electr. Electron. Eng., GE-22, 3-13. Fehler, M., L. Huang, R.S. Wu and X.B. Xie, 2005, Seismic image resolution: Numerical investigation of role of migration imaging operator, Expanded abstracts, SEG 75th Annual Meeting, 1870-1873. Fletcher, R.F., P. Fowler, P. Kitchenside, and U. Albertin, 2005, Suppressing artifacts in prestack reverse time migration: 75th Annual International Meeting, SEG, Expanded Abstracts, 2049-2052. Gelius, L.J., I., Lecomte, and H., Tabti, 2002, Analysis of the resolution function in seismic prestack depth imaging: Geophysical Prospecting, 50, 505-515. Gibson, R.L. and Tzimeas, C., 2002, Quantitative measures of image resolution for seismic survey design, Geophysics, 67, 1844-1852. Grimbergen, J. L. T., F. J. Dessing, and K. Wapenaar, 1998, Modal expansion of one–way operators in laterally varying media: Geophysics, 63, 995-1005. Han, Q., and R. S. Wu, 2005, One-way dual-domain propagators for scalar qP-wave in VTI media: Geophysics, 70, D9-D17. Hanitzsch, C., 1995, Amplitude preserving prestack Kirchhoff depth migration/inversion in laterally inhomogeneous media: Ph.D. dissertation, University of Karlsruhe. Harris, J.M., 1987, Diffraction tomography with discrete arrays of sources and receivers: IEEE Trans. Geoscience and Remote Sensing, Vol. GE-25, No. 4, 448-455. He, Y.F., and R.S. Wu, 2007, One-return boundary element method and salt internal multiples: 77th Annual International Meeting, SEG, Expanded Abstracts, 2080–2083. Hoffmann, J., 2001, Illumination, resolution, and image quality of PP- and PS-waves for survey planning: The Leading Edge, 20, 1008–1014. Hou, A., and K. Marfurt, 2005, 3D PSDW prestack depth migration: 75th Annual International Meeting, SEG, Expanded Abstracts, 1818-1821.

119

Huang, L. J., M. Fehler, and R. S. Wu, 1999a, Extended local Born Fourier migration method: Geophysics, 64, 1524-1534. Huang, L. J., M. C. Fehler, P. M. Roberts, and C.C. Burch, 1999b, Extended local Rytov Fourier migration method: Geophysics, 64, 1535-1545. Huang, L., and M. C. Fehler, 2000, Globally optimized Fourier finite-difference migration method: 70th Annual International Meeting, SEG, Expanded Abstracts, 802-805. Hubral, P., M. Tygel, and H. Zien, 1991, Three-dimensional, true-amplitude zero-offset migration: Geophysics, 56, 18-26. Jia, X.F., and R.S. Wu, 2007, Imaging steep salt flanks by super-wide angle one-way method: 77th Annual International Meeting, SEG, Expanded Abstracts, 2265-2269. Jin, S., R. S. Wu, and C. Peng, 1998, Prestack depth migration using a hybrid pseudo-screen propagator: 68th Annual International Meeting, SEG, Expanded Abstracts, 1819-1822. Jin, S., C. C. Mosher, and R. S. Wu, 2002, Offset-domain pseudoscreen prestack depth migration: Geophysics, 67, 1895-1902. Jin, S., and D., Walraven, 2003, Wave equation GSP prestack depth migration and illumination: The Leading Edge, 22, 606-610. Jin, S., M., Luo, S., Xu, and D., Walraven, 2006, Illumination amplitude correction with beamlet migration, The Leading Edge, 25, 1046-1050. Jin, S., S. Xu, and D. Walraven, 2006, One-return wave equation migration: Imaging of duplex waves: 76th Annual International Meeting, SEG, Expanded Abstracts, 2338–2342. Jo, C. H., C. Shin, and J. H. Suh, 1996, An optimal 9-point, finite-difference, frequency-space, 2D scalar wave extrapolator: Geophysics, 61, 529–537. Kiyashchenko, D., R.-E. Plessix, B. Kashtan and V. Troyan, 2005, Improved amplitude multione-way modeling method: Wave Motion, 43, 99-115. Lambar´e, G., St´ephane Opertoz, Pascal Podvin, and Philippe Thierry, 2003, 3D ray+Born migration/inversion—Part 1: Theory, Geophysics, 68, 1348-1356. Le Rousseau, J. H., and M. V. de Hoop, 2001, Scalar generalized–screen algorithms in transversely isotropic media with a vertical symmetry axis: Geophysics, 66, 1538-1550. Lecomte, I., and L.-J. Gelius, 1998, Have a look at the resolution of prestack depth migration for any model, survey and wavefields, paper presented at Annual Meeting Soc. of Explor. Geophys., New Orleans, La. Lecomte, I., Gjoystdal, H. and Drottning, A., 2003, Simulated prestack local imaging: A robust and efficient interpretation tool to control illumination, resolution, and time-lapse properties of reservoirs, 73rd Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1525-1528. Levin, S.A., 1998, Resolution in seismic imaging: Is it all a matter of perspective? Geophys., 63, 743-749. Lee, M., and S.Y. Suh, 1985, Optimization of one-way wave equations: Geophysics, 50, 1634– 1637.

120

Li, P.M., and L.G. Dong, 2006, Optimal seismic survey design based on seismic wave illumination: 68th Annual International Meeting, EAGE, Extended Abstracts, F029. Luo, M. and R.S. Wu, 2003, 3D beamlet prestack depth migration using the local cosine basis propagator: 73rd Annual International Meeting, SEG, Expanded Abstracts, 985-988. Luo, M.Q, J. Cao, X.B. Xie, and R.S. Wu, 2004, Comparison of illumination analyses using oneway and full-wave propagators: 74th Annual International Meeting, SEG, Expanded Abstracts, 67–70. Luo, M., R. S. Wu, and X. B. Xie, 2004, Beamlet migration using local cosine basis with shifting windows: 74th Annual International Meeting, SEG, Expanded Abstracts, 945-948. Luo, M.Q., R.S. Wu, and X.B. Xie, 2005, True amplitude one-way propagators implemented with localized corrections on beamlets: 75th Annual International Meeting, SEG, Expanded Abstracts, 1966-1969. Lysmer, J., and L. A. Drake, 1972, A finite-element method for seismology, in Bolt, B. A., Ed., Methods in computational physics, 11: Seismology: Surface waves and earth oscillations: Academic Press, Inc. Mao, J. and Wu, R.S., 2007, Illumination analysis using local exponential beamlets, Expanded Abstracts, SEG 77th Annual Meeting, 2235-2239. Mallat, S., 1998, A Wavelet Tour of Signal Processing: Academic Press. Mallat, S., 1999, A Wavelet Tour of Signal Processing, second edition: Academic Press. Malvar, H.S., 1992, Signal processing with lapped transforms: Artech House. Marfurt, K. J., 1984, Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations: Geophysics, 49, 533–549. Marfurt, K. J., and C. Shin, 1989, The future of iterative modeling in geophysical exploration, in E. Eisner, Ed., Handbook of geophysical exploration: I—Seismic Exploration, 21: Supercomputers in seismic exploration: Pergamon Press, 203–228. Min, D.J., C. Shin, B.D. Kwon, and S. Chung, 2000, Improved frequency-domain elastic wave modeling using weighted-averaging difference operators: Geophysics, 65, 884-895. Min, D.J. and C. Shin, 2006, Refraction tomography using a waveform-inversion backpropagation technique, Geophysics, 71, R21-R30.. Morse, P.M., and H. Feshbach, 1953, Methods of theoretical physics: McGraw-Hill Book Company. Mosher, C., and D. Foster, 2000, Common angle imaging conditions for prestack depth migration: 70th Annual International Meeting, SEG, Expanded Abstracts, 830–833. Mosher, C.C., D.J. Foster, and S. Hassanzadeh, 1997, Common angle imaging with offset plane waves: 67th Annual International Meeting, SEG, Expanded Abstracts, 1379-1382 Muerdter, D., and D. Ratcliff, 2001a, Understanding subsalt illumination through ray-trace modeling, Part 1: Simple 2D salt models: The Leading Edge, 20, 578-594.

121

Muerdter, D., M. Kelly and D. Ratcliff, 2001b, Understanding subsalt illumination through raytrace modeling, Part 2: Dipping salt bodies, salt peaks, and nonreciprocity of subsalt amplitude response: The Leading Edge, 20, 688-687. Muerdter, D., and D. Ratcliff, 2001c, Understanding subsalt illumination through ray-trace modeling, Part 3: Salt ridges and furrows, and impact of acquisition orientation: The Leading Edge, 20, 803-816. Mufti, I.R. J.A. Pita, and R.W. Huntley, 1996, Finite-difference depth migration of explorationscale 3-D seismic data: Geophysics, 61, 776-794. Mulder, W.A., and R.E. Plessix, 2004, A comparison between one-way and two-way waveequation migration: Geophysics, 69, 1491-1504. Operto, S., Gilles Lambar´ez, Pascal Podvinz, and Philippe Thierryz, 2003, 3D ray+Born migration/inversion—Part 2: Application to the SEG/EAGE overthrust experiment, Geophysics, 1357-1370. Papanicolaou, G., L. Ryzhik, and K Solna, Statistical stability in time-reversal, SIAM J. Appl. Math., 64, 1133-1155. Pratt, R.G., 1999, Seismic waveform inversion in the frequency domain, Part 1: Thoery and verification in a physical scale model: Geophysics, 64, 888-901. Pratt, R.G. and Worthington, M.H., 1988, The application of diffraction tomography to crosshole seismic data: Geophysics, 53, 1284-1294. ———, 1990, Inverse theory applied to mulitple-source cross-hole tomography I: Acoustic wave-equation method: Geophys. Prosp., 38, 287-310. Priolo, E. and Chiaruttini, C., 2003, Analytical and numerical analysis of tomographic resolution with band-limited signals, Geophysics, 68, 600-613. Prucha, B. B., and W. Symes, 1999, Angle-domain common-image gathers by wave-equation migration: 69th Annual International Meeting, SEG, Expanded Abstracts, 824–827. Qian, S. and D. Chen, 1996, Joint Time-Frequency Analysis, Methods and Applications: Prentice Hall PTR. Rickett, J.E., 2003, Illumination-based normalization for wave-equation depth migration, geophysics, 68, 1371-1379. Rickett, J., and P. Sava, 2001, Offset and angle domain common-image gathers for shot-profile migration: 71st Annual International Meeting, SEG, Expanded Abstracts, 1115–1118. Rickett, J.E. and P.C. Sava, 2002, Offset and angle-domain common image-point gathers for shot-profile migration: Geophysics, 67, 883-889. Ristow, D., and T. Rühl, 1994, Fourier finite-difference migration: Geophysics, 59, 1882–1893. Rosales, D.A., S. Fomel, B. Biondi, and P. Sava, 2008, Wave-equation angle-domain commonimage gathers for converted waves: Geophysics, 73, S17-S26. Sava, P. C., and S. Fomel, 2003, Angle-domain common-image gathers by wavefield continuation methods: Geophysics, 68, 1065–1074.

122

Schneider,W. A., and G. A.Winbow, 1999, Efficient and accurate modeling of 3-D seismic illumination: 69th Annual International Meeting, SEG, Expanded Abstracts, 633–636. Schultz, P. S., and J. F. Claerbout, 1978, Velocity estimation and downward-continuation by wavefront synthesis: Geophysics, 43, 691–714. Schuster, G.T., and Hu, J., 2000, Green’s function for migration: Continuous recording geometry, Geophysics, 65, 167-175. Sheng, J., A. Leeds, M. Buddensiek, and G.T. Schuster, 2006, Early arrival waveform tomography on near-surface refraction data, Geophysics, 71, U47-U57. Shin, C., and H. Sohn, 1998, A frequency-space 2-D scalar wave extrapolator using extended 25point finite-difference operators: Geophysics, 63, 289–296. Sjoeberg, T., Gelius, L. and Lecomte, I., 2003, 2D deconvolution of seismic image blur, 73rd Ann. Internat. Mtg.: Soc. of Expl. Geophys., 1055-1058. Sjoeberg, T.A. and Lecomte, I., 2003, 2-D deconvolution of seismic image blur, 73rd Annual International Meeting, SEG, Expanded Abstracts, 1055-1058. Snieder, R. and Lomax, A., 1996, Wavefield smoothing and the effect of rough velocity perturbations on arrival times and amplitudes, Geophys. J. Int. 125, 796-812. Spetzler, J. and Snieder, R., 2001, The effects of small-scale heterogeneity on the arrival time of waves, Geophys. J. Int. 145, 786-796. Spetzler, J., Sivaji, C., Nishizawa, O. and Fukushima, Y., 2002, A test of ray theory and scattering theory based on a laboratory experiment using ultrasonic waves and numerical simulation be finite-difference method, Geophys. J. Int. 148, 165-178. Spetzler, J. and Snieder, R., 2004, Tutorial: The Fresnel volume and transmitted waves, Geophysics, 69, 653-663. Stekl, I., and R. G. Pratt, 1998, Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators: Geophysics, 63, 1779–1794. Steinberg, B. Z., 1993, Evolution of local spectra in smoothly varying nonhomogeneous environments-local canonization and marching algorithms: Journal of the Acoustic Society of America, 93, 2566–2580. Stolt, R.H., and A.K. Benson, 1986, Seismic migration: theory and practice, handbook of geophysical exploration, vol. 5: Geophysical Press. Sun, R., and G.A. McMechanz, 2001, Scalar reverse-time depth migration of prestack elastic seismic data: Geophysics, 66, 1519-1527. Tabti, H., L. Gelius, and I. Lecomte, 2001, Discussion On “Resolution limits of migrated images” (J. Chen and G. T. Schuster, Geophysics, 64, 1046-1053). Geophysics 66, 691693. Tarantola, A., 1987, “Inverse Problem Theory, Methods for data fitting and model parameter estimation”, Elsevier, New York. Tarantola, A., 2005, Inverse problem theory and methods for model parameter estimation: SIAM (first edition, 1987).

123

Thomson, C. J., 2005, Accuracy and efficiency considerations for wide-angle wavefield extrapolators and scattering operators: Geophysical Journal International, 163, 308-323. Volker., A. W. F., G. Blacqui `ere., A. J. Berkhout‡, and L. Ongkiehong.. 2001, Comprehensive assessment of seismic acquisition geometries by focal beams—Part II: Practical aspects and examples, Geophysics, 66, 918–931 Volker, A. W. F., Blacquière, G., Berkhout, A. J., and Ongkiehong, L., 2001, Comprehensive assessment of seismic acquisition geometries by focal beams—Part II: Practical aspects and examples, Geophysics, 66, 918-931. Wang, Y. and R.S. Wu, 2002, Beamlet prestack depth migration using local cosine basis propagator: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1340-1343. Wang, Y., R. Cook, and R.S. Wu, 2003, 3D local cosine beamlet propagator: 73rd Annual International Meeting, SEG, Expanded Abstracts, 981-984. Wang, X., D. Xia, J. Li, Q. Tang and X. Jiang, 2005, Model based processing (III): pseudospace reverse time migration: 75th Annual International Meeting, SEG, Expanded abstracts, 2261-2264. Wapenaar, C. P. A., and J. L. T. Grimbergen, 1996, Reciprocity theorems for one-way wave fields: Geophysical Journal International, 127, 169–177. Wickerhauser, M.V., 1994, Adapted wavelet analysis form theory to software: A K Peters. Widess, M.B., 1982, Quantifying resolving power of seismic systems, Geophysics, 47, 11601173. Williamson, P.R., 1991, A guide to the limits of resolution imposed by scattering in ray tomography, Geophysics, 56, 202-207. Woodward, M., 1992, Wave-equation tomography: Geophysics, 57, 15-26. Wu, R. S., 1994, Wide-angle elastic wave one-way propagator in heterogeneous media and an elastic wave complex-screen method: Journal of Geophysical Research, 99, 751-766. Wu, R. S., 1996, Synthetic seismograms in heterogeneous media by one-return approximation: Pure and Applied Geophysics, 148, 155-173. Wu, R.S., F.V. Araujo and L.J. Huang, 1994, Multi-frequency backscattering tomography for constant and vertically varying backgrounds, Int. J. of Imaging Syst. & Tech., 5, 7-21. Wu, R.S., and J. Cao, 2005, WKBJ solution and transparent propagators: 67th Annual International Meeting, EAGE, Extended Abstracts, P167. Wu, R.S., and L. Chen, 2001, Beamlet migration using Gabor-Daubechies frame propagator: 63rd Annual International Meeting, EAGE, Expanded Abstracts, 74-78. Wu, R.S., and L. Chen, 2002, Mapping directional illumination and acquisition-aperture efficacy by beamlet propagators: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1352-1355. Wu, R.S. and Chen, L., 2002b, Wave propagation and imaging using Gabor-Daubechies beamlets: "Theoretical and Computational Acoustics", World Scientific, New Jersey, 661-670.

124

Wu, R.S., and L. Chen, 2003, Prestack depth migration in angle-domain using beamlet decomposition: Local image matrix and local AVA: 73rd Annual International Meeting, SEG, Expanded Abstracts, 973–976. Wu, R.S., and L. Chen, 2006, Directional illumination analysis using beamlet decomposition and propagation: Geophysics, 71, S147-S159. Wu, R.S., and L. Chen, 2006b, Target-oriented beamlet migration based on Gabor-Daubechies frame decomposition: Geophysics, 71, S37–S52. Wu, R.S., Chen, S. and Luo, M., 2004, Migration amplitude correction in angle domain using beamlet decomposition: Expanded abstracts, EAGE 66th Annual Meeting. Wu, R. S., L. Chen, and X. B. Xie, 2003, Directional illumination and acquisition dip-response, 65th Annual Conference and Exhibition, EAGE, Extended Abstracts, P147. Wu, R. S., and J. Cao, 2005, WKBJ solution and transparent propagators, 67th Annual International Meeting: EAGE, Extended Abstracts, P167. Wu, R. S., and X. Jia, 2006, Accuracy improvement for super-wide angle one-way waves by wavefront reconstruction: 76th Annual International Meeting, SEG, Expanded Abstracts, 2976-2980. Wu, R.S., M.Q. Luo, S.C. Chen, and X.B. Xie, 2004, Acquisition aperture correction in angledomain and true-amplitude imaging for wave equation migration: 74th Annual International Meeting, SEG, Expanded Abstracts, 937-940. Wu, R.S., and J. Mao, 2007, Beamlet migration using local harmonic bases: 77th Annual International Meeting, SEG, Expanded Abstracts, 2230-2234. Wu, R.S. and Toksöz, M.N., 1987, Diffraction tomography and multi-source holography applied to seismic imaging, Geophysics, 52, 11-25. Wu, R.S., X.B. Xie, M. Fehler and L.J. Huang, 2006, Resolution analysis of seismic imaging: 68th Annual International Meeting, EAGE, Expanded Abstracts, G048. Wu, R.S., Y. Wang and J.H. Gao, 2000, Beamlet migration based on local perturbation theory: 70th Annual Meeting, SEG, Expanded Abstracts, 1008-1011. Wu, R.S., Y. Wang and M.Q. Luo, 2003, Local-cosine beamlet migration for 3D complex structures: 8th International Congress of the Brazilian Geophysical Society, 14-18. Wu, R.S., Y. Wang and M.Q. Luo, 2008, Beamlet migration using local cosine basis: Geophysics, accepted. Wu, W.J., L.R. Line, and H.X. Lu, 1996. Analysis of higher-order, finite-difference schemes in 3-D reverse-time migration: Geophysics, 61, 845-856. Wu, R.S. ,Wu, Y.L., Liu, C.R., Peng, C.Y., and M.J. Wang, 1977, Multi-frequency synthetic detecting holography, Acta Geophysica Sinica, 20, 33-38. Xie, X. B., Z., Ge and T., Lay, 2005, Investigating explosion source energy partitioning and Lgwave excitation using a finite-Difference plus slowness analysis method, Bull. Seism. Soc. Am. 95, 2412-2427.

125

Xie, X. B., S.W. Jin, and R. S. Wu, 2003, Three-dimensional illumination analysis using waveequation based propagator: 73rd Annual International Meeting, SEG, Expanded Abstracts, 989–992. Xie, X.B., S.W. Jin and R.S. Wu, 2004, Wave equation based illumination analysis: 74th Annual International Meeting, SEG, Expanded Abstracts, 933-936. Xie, X.B., S.W. Jin and R.S. Wu, 2006, Wave-equation-based seismic illumination analysis: Geophysics, 71, S169-S177. Xie, X. B., C. C. Mosher, and R. S. Wu, 2000, The application of wide angle screen propagator to 2D and 3D depth migrations: 70th Annual International Meeting, SEG, Expanded Abstracts, 878–881. Xie, X.B. and R.S. Wu, 1998, Improve the wide angle accuracy of screen method under large contrast: 68th Annual International Meeting, SEG, Expanded Abstracts, 1811-1814. Xie, X. B., and R. S. Wu, 2005, Multicomponent prestack depth migration using elastic screen method: Geophysics, 70, S30-S37. Xie, X.B., and R.S. Wu, 2002, Extracting angle domain information from migrated wave field: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1360-1363. Xie, X.B., R.S. Wu, M. Fehler and L. Huang, 2005, HSeismic resolution and illumination: A wave-equation-based analysisH, 75th Annual International Meeting, SEG, Expanded Abstracts, 1862-1865. Xie, X.B., and R.S., Wu, 2006, A depth migration method based on the full-wave reverse-time calculation and local one-way propagation, 76th Annual International Meeting, SEG, Expanded Abstracts, 2333-2337. Xie, X.B., Wu, R.S., M.Fehler and L. Huang, 2005, Seismic resolution and illumination: A wave-equation-based analysis: 75th Annual International Meeting, SEG, Expanded abstracts, 1862-1865. Xie, X.B., and H., Yang, 2008, A full-wave equation based seismic illumination analysis method, 70th Conference and Technical Exhibition, EAGE, Expanded abstracts, accepted. Xu, S., H. Chauris, G. Lambaré, and M. Noble, 2001, Common-angle migration: A strategy for imaging complex media: Geophysics, 66, 1877-1894. Xu, S., and S. Jin, 2006, Wave equation migration of turning waves: 76th Annual International Meeting, SEG, Expanded Abstracts, 2328–2332. Xu, T., McMechan, G.A. and Sun, R., 1995, 3-D prestack full-wavefield inversion: Geophysics, 60, 1805-1818. Yang, H., X.B. Xie, M.Q. Luo, and S.W. Jin, 2008, Target oriented full-wave equation based illumination analysis: 78th Annual International Meeting, SEG, Expanded Abstracts, submitted. Yoon, K., C. Shin, S. Suh, L. Lines, S.H., Kigam, 2003, 3D reverse-time migration using the acoustic wave equation: An experience with the SEG/EAGE data set: The Leading Edge, 22, 38-41.

126

Yoon, K., K.J. Marfurt, and W. Starr, 2004, Challenges in reverse-time migration: 74th Annual International Meeting, SEG, Expanded Abstracts, 1057-1060. Yu, J., and Schuster, G.T., 2003, 3-D prestack migration decomvolution: 73rd Annual International Meeting, SEG, Expanded Abstracts, 1651-1654. Zhou, C., W. Cai, Y. Luo, G.T. Schuster, and S. Hassanzadeh, 1995, Acoustic wave-equation traveltime and waveform inversion of crosshole seismic data: Geophysics, 60, 765-773. Zhou, C., G.T. Schuster, S. Hassanzadeh, and J.M. Harris, 1997, Elastic wave-equation traveltime and waveform inversion of crosswell seismic data: Geophysics, 62, 853-868. Zhang, G.Q., 1993, System of coupled equations for upgoing and downgoing waves: Acta Mathematicae Applicatae Sinica, 16, 251–263. Zhang, Y., G. Zhang, and N. Bleistein, 2003, True amplitude wave equation migration arising from true amplitude one-way wave equations: Inverse Problems, 19, 1113-1138. Zhang, Y., G. Zhang, and N. Bleistein, 2005, Theory of true-amplitude one-way wave equations and true-amplitude common-shot migration: Geophysics, 70, E1–E10. Zhang, Y., S. Xu, and G.Q. Zhang, 2006, Imaging complex salt bodies with turning-wave oneway wave equation: 76th Annual International Meeting, SEG, Expanded Abstracts, 2323–2326.

127