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
RuShan Wu and XiaoBi 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. DEFG0204ER15530
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. AngularSpectral 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 WAVEEQUATIONBASED 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: TRUEREFLECTION IMAGING .................................................................................................................................... 35 4. OPTIMIZATIONS IN TRUEAMPLITUDE ONEWAY PROPAGATORS..................... 36 4.1. Implementations of the trueamplitude oneway wave equation ................................... 37 4.2. Optimizations for different terms in the trueamplitude oneway propagators ............. 39 4.3. Numerical results ........................................................................................................... 41 4.4. Conclusions .................................................................................................................... 45 5. SUPERWIDEANGLE ONEWAY 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 FULLWAVE 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. LOCALANGLE DOMAIN ILLUMINATION IN FREQUENCY DOMAIN FOR FULLWAVE PROPAGATORS......................................................................................................... 68 7.1. Wavefield Decomposition for OneWay Propagators ................................................... 69 7.2. Wavefield Decomposition for FullWave Propagators.................................................. 70 7.3. LocalAngle Domain Illumination Results .................................................................... 72
2
7.4. Discussion About The Influence of Multiples on The Illumination Strength ............... 76 7.5. Conclusions .................................................................................................................... 78 Appendix A GaborDaubechies 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: GaborDaubechies Frame (GDF) Beanlet Decomposition ............................. 97 Appendix B: Local Exponential Frame (Lef) Beamlet Decomposition ............................... 97 9. A DEPTH MIGRATION METHOD BASED ON THE FULLWAVE REVERSETIME CALCULATION AND LOCAL ONEWAY 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 AngleDomain 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, 19461949. 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, 24992503. Cao, J., and R.S. Wu, 2008, Amplitude compensation for oneway wave propagators in inhomogeneous media and its application to seismic imaging, Communications in Computational physics, 3, 203221. 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, 18701873. He, Y.F., and Wu, R.S., 2007, Onereturn boundary element method and salt internal multiples, 77th Annual International Meeting, SEG, Expanded abstracts, 20802084. Jia, X.F., and Wu, R.S., 2007, Imaging steep salt flanks by superwide angle oneway method, 77th Annual International Meeting, SEG, Expanded Abstracts, 22652269. Luo, M., R.S. Wu, and X.B. Xie, 2005, True amplitude oneway propagators implemented with localized corrections on beamlets, 75th Annual International Meeting, SEG, Expanded Abstracts, 19661969. Mao, J., and Wu, R.S., 2007, Illumination analysis using local exponential beamlets, Expanded abstracts, 77th Annual International Meeting, SEG, Expanded Abstracts, 22352239. Wu, R.S., and Jia, X.F., 2006, Accuracy improvement for superwide angle oneway waves by wavefront reconstruction, 76th Annual International Meeting, SEG, Expanded abstracts, 29762980. Wu, R.S., X.B. Xie, and X.Y. Wu, 2006, oneway and onereturn approximations (dewolf approximation) for fast elastic wave modeling in complex media, in: Advances in Wave Propagation in Heterogeneous Earth, Elsevier, 265322. Wu, R.S., X.Y. Wu, and X.B. Xie, 2006, Simulation of highfrequency wave propagation in complex crustal waveguides using generalized screen propagators, in: Advances in Wave Propagation in Heterogeneous Earth, Elsevier, 323363. 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, 27282732. Wu, R.S., and Mao, J., 2007, Beamlet migration using local harmonic bases, 77th Annual International Meeting, SEG, Expanded abstracts, 22302234. Wu, R.S., Wu, B.Y. and Geng, Y. 2008, Seismic wave propagation and imaging using timespace wavelets, Expanded abstracts, SEG 78th Annual Meeting, 29832987. Xie, X.B., R.S.,Wu, M. Fehler and L. Huang, 2005, Seismic resolution and illumination: A waveequationbased analysis, 75th Annual International Meeting, SEG, Expanded Abstracts, 18621865. Xie, X.B., S.W. Jin and R.S. Wu, 2006, Wave equation based seismic illumination analysis, Geophysics, 71, S169177. Xie, X.B. and R.S. Wu, 2006, A depth migration method based on the fullwave reversetime calculation and local oneway propagation, 76th Annual International Meeting, SEG, Expanded Abstracts, 23332337. Xie, X.B. and H. Yang, 2007, A migration velocity updating method based on the shot index common image gather and finitefrequency sensitivity kernel, 77th Annual International Meeting, SEG, Expanded Abstracts, 27672771. Xie, X.B., and H. Yang, 2008, A fullwave equation based seismic illumination analysis method, 70th Conference & Technical Exhibition, EAGE, Expanded Abstracts, P284. Xie, X.B., and H. Yang, 2008, The finitefrequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis, Geophysics, 73, S241S249. Zhu, X., and R. Wu, 2008, Imaging diffraction points using the local image matrix in prestack migration, 78th Annual International Meeting, SEG, Expanded Abstracts, 23612365.
5
ABSTRACT Our proposed work on high resolution/high fidelity seismic imaging focused on three general areas: (1) development of new, more efficient, waveequationbased propagators and imaging conditions, (2) developments towards amplitudepreserving 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 truereflection imaging. Truereflection 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 reflectionangle. 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 acquisitionaperture effect. We have made significant progress in both categories. We studied the effects of different terms in the trueamplitude 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 oneway propagators and fullwave propagators. We developed the fast acquisitionaperture correction method in the local angledomain, which is an important element in the truereflection imaging. Other developments include the superwide angle oneway propagator and special fullwave reversetime migration method. Finally, we studied the theoretical basis of truereflection 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 waveequation 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 waveequation 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, waveequationbased propagators and imaging conditions, (2) developments towards amplitudepreserving 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 angledomain decomposition to the resolution operator and derived the angularspectral 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 acquisitionsystem 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”. Truereflection imaging is an advanced imaging technology which aims at keeping the image amplitude proportional to the reflection strength of the local interface. Truereflection image can present the image as the total strength or the reflections as a function of reflectionangle. 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 acquisitionaperture effect. In Chapter 4 we discuss the waveequation based trueamplitude propagator. We discussed the different terms in the trueamplitude oneway 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 superwide angle oneway 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 oneway wave propagator, there is a trend in the industry to use the fullwave propagator for seismic imaging, even with significantly increased cost. To apply the concepts of truereflection imaging to the fullwave imaging methods, such as
7
the finite difference reversetime migration, we need to analyze the acquisitionaperture effects using directional illumination analysis. In Chapters 6 and 7 we summarize our results on directional illumination analyses using timespace 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 fullwave migration such as reversetime migration. In Chapter 8 we present the results of our study on fast acquisitionaperture correction method in truereflection 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 acquisitionaperture correction. This paved the way for practical application of truereflection imaging in the exploration industry. In Chapter 9 a depth migration method of fullwave reversetime migration which uses local oneway 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 truereflection 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 raytracing 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 highfrequency 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 wellposed problem, we should have RI (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 offdiagonal 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 zerotime 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 (trueamplitude 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 spacedomain 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 crosscorrelation 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 dipangle 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 zerooffset (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 dk 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 explodingreflector 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 zerooffset imaging process or imaging with explodingreflector 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. fullwave FD propagators, oneway FD propagators, oneway dualdomain 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. AngularSpectral Representation of Point Spreading Function (PSF) Angularspectral representation of resolution of PSF is more intuitive. We can see directly the information coverage in local angledomain. 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 planewave 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 highfrequency 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 zaxis; 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 angledomain (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 angledomain 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 saltbody 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 zerooffset and exploding reflectors modeling respectively: (K k ; x 0 ) 2 dk 2 dxs GI * ( , x 0 , k s ; xs )GM ( , x 0 ; xs ) As
(K k s ; x0 ) 2 dk 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, traveltime variation, amplitude changes of the wavefield (forward scattering), but only sharp boundaries produce significant backscattering (reflection) and largeangle 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 oneway 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 spacedomain 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 singlefrequency. 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 angledomain 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 hf asymptotic sphere (a circle in the 2D case). The resolution is only for the dipdirection (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 singledip 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 dipangles 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 singlefrequency
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) reflectionangle 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) reflectionangle 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 ) dk 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 dipresponse (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 oneway 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) oneway wave propagator. For the acquisition resolution, only the oneway 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 oneway propagator.
22
Figure 2.9 shows the singlefrequency acquisition resolution of the acquisition system for the case of SEG/EAGE salt model. The propagator is the oneway 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 multifrequency 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 multifrequency acquisition (f = 624 Hz) for a reflection element with dipangle 45 degrees for a common scatteringangle 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 offshore acquisition using streamer receiver lines with shot in the righthand side of the streamer.
Figure 2.9: Singlefrequency 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: Multifrequency 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 multifrequency acquisition (f = 624 Hz) for a reflection element with dipangle 45 degrees for a common scatteringangle 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 oneway migration and rayKirchhoff 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 oneway waveequation 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 rayKirchhoff 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 waveequation imaging; on the bottom right is the PDF from the rayKirchhoff imaging.
Figure 2.14: Lateral resolutions at depth 2 km (see Figure 2.7) of waveequation imaging and rayKirchhoff imaging in a random medium.
27
3. SEISMIC RESOLUTION AND ILLUMINATION: A WAVEEQUATIONBASED 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 raybased highfrequency 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 waveequationbased 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 waveequationbased 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 oneway 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: 2D 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 receiverside 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 40degree 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 waveequation propagator, a group of numerical examples are calculated using the 2D SEG/EAGE salt model and its acquisition configuration. The wavefield is extrapolated using a wideangle oneway 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 subsalt region, the resolution function is weak and distorted, resulting in poorquality 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 40degree 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 waveequationbased 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: TRUEREFLECTION IMAGING
35
4. OPTIMIZATIONS IN TRUEAMPLITUDE ONEWAY PROPAGATORS In the traditional oneway 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 raytheoretical amplitudes (Zhang, et al., 2003). WKBJ amplitude was then introduced into the original oneway 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 fullwave operator into oneway wave operators in heterogeneous media (Zhang, 1993). With an extra amplitude term introduced to the conventional oneway wave equations, the firstorder transport equation of the oneway wave equation is the same as that from the fullwave equation in the sense of highfrequency approximation. In this sense, they are called “trueamplitude” oneway 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 fluxnormalized propagator for general heterogeneous media (Wapenaar & Grimbergen, 1996). In the traditional oneway propagators, the square root operator is approximated by a series of rational functions. One approach to achieve better accuracy for the large propagationangle waves is to use higherorder 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 oneway ‘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 leastsquares optimization scheme. Ristow & Rühl (1994) proposed a locally optimized scheme for the Fourier FiniteDifference (FFD) propagator. Huang & Fehler (2000) proposed a globally optimized scheme for the WideAngle Screen (WAS) propagator or GSP (Generalized Screen Propagator). Considering the advantage of dualdomain 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 oneway extrapolator here. In this paper, three implementations for the amplitude correction in the trueamplitude oneway 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 propagationangle 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 trueamplitude oneway wave equation The trueamplitude oneway 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 squareroot operator for the traditional oneway 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 oneway 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 oneway wave propagator (4.4), firstorder 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 firstorder 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 dualdomain 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).
Onestep 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 oneway wave eq. (4.4) and the wideangle FD correction in WAS and FFD propagator. Their computation costs should be similar. Twostep WKBJ correction in dual (space and wavenumber) domain The advantages of above onestep FD implementation for WKBJ correction are that no reference velocity is required which makes this implementation method compatible with any traditional oneway extrapolator; and the FD correction method in space domain should be stable. However, considering the disadvantage of pure FD extrapolators compared with the dualdomain extrapolators, an alternative dual domain twostep 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 onestep FD correction (4.9), just substituting with 0. This dual domain twostep WKBJ correction can easily be implemented in the traditional dual domain oneway propagators (e.g., WAS, FFD). 0 z
2
2 x
2 y
Targetoriented 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 timeconsuming if we are interested in the amplitude behavior of the pressure only in a very deep target, e.g. subsalt 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 oneway 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 energyflux 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 targetoriented 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 energyflux 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 firstorder 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 trueamplitude oneway propagators
39
Different oneway 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 propagationangle waves if the optimized coefficients for (4.6) are applied to (4.8) (Figure 4.1). Compared with traditional oneway 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 oneway (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 firstorder 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 firstorder Padé approximation to , eq. (4.6), we use the expansion coefficients derived from the leastsquares optimization scheme by Lee & Suh (1985). After the optimization, the relative errors for above 3 approximated operators are much smaller for the larger propagationangle 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 firstorder 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: firstorder Padé approximation to the square root
40
operator (4.6) (red lines), WAS approximation to the square root operator (4.8) (green lines), and firstorder 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 oneway propagator. Then above discussed three optimized implementations, onestep FD method, twostep dual domain method and flux transparent method, are compared in a smoothly varying heterogeneous model. Finally the WKBJcorrected amplitude is investigated in the smoothed 2D SEG/EAGE salt model. Phase distortion caused amplitude error In the traditional oneway propagator, much less attention is paid to the amplitude than to the phase. Before the amplitude WKBJ correction is applied to the traditional oneway propagator, we first investigate the amplitude behavior related to the phase errors in the traditional oneway 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 phasescreen and WAS propagator compared with the true values. For phasescreen propagator (Figure 4.3a, b), the phase for 0º propagationangle 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 propagationangle bigger than 0º, which distorts the energy flux of the 0º propagationangle wave. With the wideangle phase more accurate WAS propagator, the phase of 60º propagationangle 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 oneway 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 propagationangle waves, optimization should also be applied to the traditional oneway 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 phasescreen 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 trueamplitude 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/xs1/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 onestep FD correction in space domain the amplitude is corrected very well up to 60º (Figure 4.5b). With the twostep method the amplitude is corrected very well up to 75º (Figure 4.5c), which is better than the onestep 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 WKBJcorrected 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 subsalt area (green line in Figure 4.6a). The amplitudes from onestep 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 WKBJcorrected amplitude since only the pressure at the receivers’ depth level needs to be recovered from the downwardcontinued flux field.
42
Figure 4.3: Phase and amplitude behavior for the phasescreen and WAS propagator with a constant velocity contrast m=1.2. (a), (b) are for phasescreen 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 phasescreen 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.121/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.121/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.121/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.121/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.121/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/xs1/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 fullwave FD method, and the dashed lines are results from oneway 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 twostep 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 subsalt area (green line); (b) picked maximum amplitude with/out the WKBJ correction along the receivers shown in (a).
4.4. Conclusions Three implementations (onestep FD method, twostep dual domain method and targetoriented flux transparent method) for the amplitude correction in TAOP are investigated. In the traditional oneway 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 propagationangle waves, optimization schemes are applied to all the approximated terms in TAOP including the traditional oneway (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 oneway extrapolators: twostep method is more accurate than onestep method but onestep method is compatible with any traditional oneway extrapolator; the flux transparent method is less accurate than onestep 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 onestep method; the advantage of the flux transparent method is that targetoriented 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 fullwave FD method, and there is almost no extra cost for the flux transparent method to obtain this WKBJcorrected amplitude for that target.
46
5. SUPERWIDEANGLE ONEWAY MIGRATION: APPLICATION TO 2004 BP 2D BENCHMARK MODEL 5.1. Background The oneway wave equation simulates the exact oneway wavefield only in very limited propagation directions. Quite a few methods have been proposed to improve the wideangle accuracy of oneway 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 wideangle waves for strong contrast media is still a serious problem and put a practical limitation on applying these methods to steep reflector imaging. The superwideangle oneway method (Wu and Jia, 2006) has been developed to improve the regular oneway propagators by a wavefront reconstruction scheme which combines and interpolates the weighted orthogonally propagated oneway 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 oneway wave equations in z and xdirections 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 oneway wavefields propagated downward in z
(5.1)
direction and horizontally in xdirection, respectively; PD and PH are the corresponding oneway 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 righthand side of eq. (5.2). Now considering the oneway propagations in zdirection having n steps forward and in xdirection having m steps forward, we have uD x, z0 nz PDn uD x, z0 , uH x0 mx, z PHm uH x0 , z .
47
(5.3)
We reconstruct the accurate wavefield by taking a weighted average of the two oneway wavefields. Wavefront reconstruction is a stepwise marching process, which means the wavefront reconstructed at the current step will be used for further oneway propagation. Note that the oneway propagation is exact in the forward direction and very accurate for smallangle waves. Therefore, at any location we put heavy weight on the oneway solution which propagates smallangle waves at that point and light weight on the other oneway solution which propagates largeangle 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 xaxis; 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 oneway 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 xdirections, respectively. This corresponds to a simple scheme of weighted average without updating the wavefront interactively. The interpolation of the two oneway wavefields is separated completely from the wave propagations. As a result, the propagations remain the same as the regular oneway 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 oneway wavefields. We can use either horizontal screens (i.e. downward oneway wave propagation) or vertical screens (i.e. horizontal oneway 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 oneway wave propagations.
5.3. Experiments on 2004 BP 2D Benchmark Dataset The weighting summation can be applied to most of the oneway propagators. In this work, we test the scheme using the wideangle 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 superwideangle oneway migration will be demonstrated on the 2004 BP 2D dataset (Billette and BrandsbergDahl, 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 oneway migration methods.
Figure 5.1: velocity model of 2004 BP 2D benchmark dataset.
Figure 5.2: modeling results of the superwideangle 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 superwideangle 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 superwideangle GSP method. Figure 5.4 shows the singleshot 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 oneway method for migration. However, in the imaging results of the superwideangle 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 superwideangle 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 oneway pass in addition to the regular superwideangle 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 oneway wavefields. Wavefield gradients can still be used to determine the weights in this case. When the superwideangle oneway 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 superwideangle oneway 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 oneway 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 oneway propagation methods.
(a)
(b) Figure 5.4: comparison of the singleshot images obtained by the regular GSP method and the superwideangle GSP method with the shot location at x=28km. (a) regular downward GSP method; (b) superwideangle GSP method.
51
(a)
(b) Figure 5.5: comparison of the stacked images obtained by the regular GSP method and the superwideangle GSP method with the shots from x=25km to x=30km. (a) regular downward GSP method; (b) superwideangle GSP method.
52
(a)
(b) Figure 5.6: Comparison of the singleshot images obtained by the regular LCB method and the superwideangle LCB method with the shot location at x=28km. (a) regular downward LCB method; (b) superwideangle LCB method.
53
(a)
(b) Figure 5.7: comparison of the stacked images obtained by the regular LCB method and the superwideangle LCB method with the shots from x=25km to x=30km. (a) regular downward LCB method; (b) superwideangle LCB method.
54
Figure 5.8: image obtained by superwideangle GSP plus downward GSP method.
55
6. TARGET ORIENTED FULLWAVE 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 incidentscattering waves from different directions. The raybased 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 highfrequency asymptotic approximation causes intrinsic difficulties for this method being used in complex overburden structures. Recently, the oneway 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 oneway 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 wideangle signals such as turning waves. In the last few years, reversetime 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 reversetime migration. To overcome the angle limitation of the oneway propagator, Xie and Yang (2008) proposed a method which uses the fullwave finitedifference method as the propagator and uses a timedomain 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 splitstep local planwave 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 fullwave propagator. The new method does not 56
have angle limitations usually encountered by a oneway approach. In addition, it accurately solves the energy transmission in a complex velocity model. This method is particularly useful in providing illumination analysis for reversetime 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 incidencescattering beams forms a basic illumination measurement (see Figure 6.1). Note, due to multipathing, 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 sourcereceiver pairs A r , s , g A r , s , g ; rs , rg . (6.5) rs rg
where Ar, 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 overhang 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 fullwave illumination method (the entire outer square) is much larger than that can be handled by a oneway propagator based illumination method (the inner square). This makes the fullwave 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 sourcereceiver pair, (b) sourceside and (c) receiverside 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 sourcereceiver 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 sourcereceiver 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 sourcereceiver 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 oneway 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.7ad are ADRs for 135, 90, 45, and 0 degrees in a v z model and generated by a pair of sourcereceiver. 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 20Hz 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 sourcereceiver pair.
Figure 6.7, ADR maps generated by a sourcereceiver 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 sourcereceiver pair in a twolayer model: (a) 45º ADR for downdown waves; (b) 120º ADR for upup waves; and (c) 120º ADR for downup waves.
Figure 6.9: Comparison of ADR maps from different kind of waves by one sourcereceiver 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 downdown waves; and (c) 120º ADR for downup waves. (d) 120º ADR for full waves. The FD calculated Green’s functions compose of direct waves, reflections, turning waves and multiples, etc. The fullwave illumination analysis includes contributions from all these waves, while a oneway propagator usually cannot handle upgoing waves. Shown in Figure 6.8 are ADRs for a twolayer 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 “downdown” waves. Figure 6.8b is the 120o ADR generated by “upup” reflections. Figure 6.8c is the 120o ADR obtained from the direct source wave and reflected receiver wave, i.e. the “downup” waves. The later two ADRs cannot be properly handled by a oneway 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 9bd 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.9bd. 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.9ad. Illumination as a function of target reflection angle: For a locally planar target with a normal vector pointing to direction n , the targetoriented 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 reversetime image for a salt dome model. Figure 6.10a shows the reversetime 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 reversetime image and target illumination for the salt dome model with (a) reversetime 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.11cd 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 BrandsbergDahl, 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 oneway wave equation based method, while the reverse time migration can image this structure satisfactorily. The illumination method based on the fullwave 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 fullwave equation based illumination analysis method which has the target oriented capability. The fullwave 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. LOCALANGLE DOMAIN ILLUMINATION IN FREQUENCY DOMAIN FOR FULLWAVE 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 raybased 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 raybased method is convenient and can provide both intensity and directional information carried in the wavefield. However, the highfrequency 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 overprecise 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 frequencydependent illumination, we need a wavetheory based method. Oneway wave equation based propagators are widely used in migration and illumination analysis. Although they neglect multiples, they properly handle multiforwardscattering 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 localizedangle 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), GaborDaubechies 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 oneway 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 oneway wave propagators have been studied by different approaches for a long time, e.g., the “trueamplitude” oneway wave equation method (Zhang, 1993; Zhang et al., 2003, 2005), multioneway method (Kiyashchenko et al., 2005). However, even with these corrections, the oneway propagator still cannot give accurate wavefield amplitude in complex models with sharp contrast. The numerical implementations based on the oneway wave equation with zaxis as the preferred propagation direction always have inherent limitation in wideangle accuracy. Some oneway propagation based methods can model the 68
backscattered 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 fullwave propagators is very desirable as a complementary method. The fullwave propagator, e.g., finitedifference (FD) and finiteelement 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 timedomain methods and the fact that the illumination is frequencydependent, we propose to analyze the illumination in the frequency domain. We can use frequency domain fullwave 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 oneway wave propagator and fullwave propagator (where single frequency full waves are extracted from the results by timedomain FD). However, artificial interference patterns exist in the illumination results from the fullwave 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 fullwave 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 splitstep method using 1D/2D decomposition to obtain LAD illumination in the frequency domain for 2D/3D fullwave propagators. Firstly we summarize the LSS and beamlet decomposition techniques for oneway propagators; based on these we give the wavefield decomposition methods for fullwave 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 OneWay 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 oneway propagators but not for fullwave propagators. Firstly we summarize these two decomposition methods for oneway propagators. Local slant stack for oneway 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 oneway 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 FullWave 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 oneway 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 upgoing 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 fullwave 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. Splitstep decomposition for fullwave propagators 2D decomposition can directly obtain the wavefield for all directions; however, it costs much more computation than 1D decomposition. Here we propose a splitstep decomposition method to obtain the LAD wavefield for fullwave propagators using more efficient 1D decomposition technique. Firstly, we decompose the fullwave 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 downgoing 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 downgoing 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. LocalAngle 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 twolayer model. Downgoing waves illumination in 2D SEG/EAGE salt model The whole acquisition system of this model consists of 325 shots with 176 leftside 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 splitstep 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 splitstep 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 (ac) and 2D LSS (df) 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 (ac) and 2D LSS (df) 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 oneway 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 fullwave 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 fullwave 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 lowerright 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 twolayer 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 splitstep 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 multiarrivals. Figure 7.10 shows the DI maps for the most energetic waves. Comparison with the DIs for the full downgoing waves (Figure 7.4ac) shows that the other arrivals (multiarrivals and multiples) provide extra illumination to the subsurface. The migration methods based on oneway propagators utilize not only the first arrival but also the multiarrivals. Therefore the multiarrivals should be included in the illumination analysis for survey design and truereflection imaging correction. However, multiples especially the internal multiples should be eliminated since most of the migration methods based on the oneway 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 mostenergetic arrival from other arrivals by timedomain windowing. However, we cannot separate multiples from multiarrivals 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 oneway and onereturn 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 subsalt 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 subsalt; however, the 77
multiples illuminate the whole subsalt more evenly and provide extra illumination to the shadow in the illumination by primaries (right part in subsalt).
7.5. Conclusions We propose an efficient splitstep local plane wave decomposition method to obtain the localangle domain (LAD) illumination in the frequency domain for fullwave 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 truereflection 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 subsalt 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 frequencydependent illumination and the frequencydomain fullwave forward modeling is usually efficient and storagesaving to provide the wavefield for given frequencies compared with the timedomain 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 fullwave modeling. We conclude that the proposed splitstep local plane wave decomposition method for fullwave propagators provides a flexible, efficient and accurate tool for the LAD wavetheory based illumination analysis in complex 2D and 3D models.
Appendix A GaborDaubechies 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 )eim x , (A7.1)
where m m ( is the wavenumber sampling interval), and g x is a Gaussian window function, x2 (A7.2) g ( x ) ( s 2 ) 1/ 4 exp 2 , 2s s2
R N 2 , 2
(A7.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 (A7.1) into (7.3), we can obtain the partially reconstructed local plane waves, (A7.4) u z ( x, m , ) eim x g ( x xn )uˆ z ( xn , m , ) , n
with the decomposition coefficients * uˆ z ( xn , m , ) uz ( x, ), bmn ( x ) uz ( x, ) bmn ( x ) dx ,
(A7.5)
where , stands for inner product, “*” stands for complex conjugate, and bmn ( x ) are the dual GD frame atoms, bmn ( x ) g ( x xn )eim x , (A7.6) with g ( x ) as the dualwindow function of g ( x ) . The dualwindow function can be calculated by pseudoinversion of the original window function (Qian and Chen, 1996; Mallat, 1998; Wu and Chen, 2001). From equation (A7.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 (A7.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 angledomain, which makes raybased methods more convenient than waveequation based methods since the angle information is inherently embedded in raybased methods. The theory and method of truereflection imaging has been developed based on highfrequency 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; BrandsbergDahl et al., 2003). However, the results may contain large errors in complex areas due to the highfrequency approximation and singularity problems therein. For waveequation 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 angledependent reflection amplitudes. This demonstrated the significance of acquisition aperture correction in truereflection imaging. For waveequation based migration methods, the subsurface angledomain commonimage gathers (CIGs) can be obtained either by decomposing the propagated wavefield before imaging (“dataspace 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 (“imagespace 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 imagespace method, the subsurface offsetdomain CIGs are firstly computed and then transformed to the subsurface angledomain CIGs. For sourcereceiver migration, the subsurface offsetdomain CIGs are naturally available. For shot profile migration, the subsurface offsetdomain 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 horizontaloffset CIGs degenerate, and their focusing around zero offset blurs (Biondi and Symes, 2004). This limitation of subsurface horizontaloffset CIGs can be sidestepped by computing subsurface offsetdomain CIGs along the vertical subsurfaceoffset axis (Biondi and Shan, 2002; Biondi and Symes, 2004). Sava and Fomel (2003) presented a simple method for transforming subsurface horizontaloffset CIGs into subsurface angledomain CIGs by a slant stack transformation (Schultz and Claerbout, 1978) applied independently to each subsurface horizontaloffset CIG. Biondi and Symes (2004) proposed a simple and effective method to create a subsurface angledomain CIG cube by combining a subsurface horizontaloffset CIG cube with a subsurface verticaloffset CIG cube. In dataspace methods, angle gathers are obtained by decomposing the propagated wavefield before imaging. The angledomain images can be obtained for shotprofile migration by de Bruin et al. (1990) using global Fourier transform in 1D media, and for sourcereceiver 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 shotprofile, sourcereceiver, or other propagators); 2. with the decomposed local plane waves we can obtain not only the angledomain image but also many other very useful angledomain 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 dipangle domain can be five times of that to obtain the spacedomain 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 angledomain imaging (Wu and Chen, 2002; 2003; Chen et al., 2006) and other anglerelated 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 frequencyspace 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 spacedomain imaging condition (Claerbout, 1971) for a single frequency is in the form of crosscorrelation of the incident sourceside wave uS ( x, ) and the backpropagated receiverside 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 spacelocal 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 truereflection 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 timeconsuming. In the original LAD implementation with beamlet decomposition, for the local plane wave with
Figure 8.1: Localangle definition in a scattering experiment for a local planar reflector. θs and θg are the local source and receiving angles, respectively. θn is reflectornormal 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 backpropagated wavefields for every shot gather; therefore, it is very timeconsuming. 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 spacelocal 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 backpropagated 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 shotsummed 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 multifrequency 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 reflectionangle 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 leftside 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 dipangle 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 dipangle 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 stripelike 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 dipangle 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 spacedomain 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 spacedomain 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 dipangle 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 stripelike 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 stripelike artifacts are removed in the amplitude correction factor for dipangle +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 dipangle +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 stripelike 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 dipangle +30º by the GDF decomposition (Figure 8.14a) now displays similar stripelike 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 stripelike 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 (A2) and (A3)) 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 dipangle +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 localangle 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 GaborDaubechies 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 spacedomain imaging. This makes it possible for the waveequation based space localangle 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 stripelike 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: GaborDaubechies 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 )eim x , (A8.1)
where m m ( is the wavenumber sampling interval), and g x is a Gaussian window function,
x2 (A8.2) g ( x ) ( s 2 ) 1/4 exp 2 , 2s R N 2 s2 , (A8.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 (A1) into (8.2), we can obtain the partially reconstructed local plane waves, u( x, m , ) eim x g ( x xn )uˆ ( xn , m , ) , (A8.4) X
X
n
with the decomposition coefficients * uˆ( xn , m , ) u( x, ), bmn ( x ) u( x, ) bmn ( x ) dx , where
,
(A8.5)
stands for inner product, “*” stands for complex conjugate, and bmn ( x ) are the dual
GD frame atoms,
bmn ( x ) g ( x xn )eim x , (A8.6) with g ( x ) as the dualwindow function of g ( x ) . The dualwindow function can be calculated by pseudoinversion of the original window function (Qian and Chen, 1996; Mallat, 1999; Wu and Chen, 2001). From equation (A8.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 , ) . (A8.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
(B8.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) , leftpropagating 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 im ( x xn ) bmn ( x ) ibmn ( x) Ln
,
(B8.2)
(c) 2 Bn ( x ) cos m ( x xn ) bmn ( x ) Ln , (B8.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 (B8.2) into (8.2), we can obtain the partially reconstructed right and leftpropagating local plane waves u ( ) ( x, m , ) and u ( ) ( x, m , ) , () 2 Bn ( x )uˆ ( ) ( xn , m , ) exp( im xn ) u ( x, m , ) exp( im x ) Ln n , 2 u ( ) ( x, , ) exp( i x ) () n L Bn ( x)uˆ ( xn , m , ) exp( im xn ) m m n with the decomposition coefficients () uˆ ( ) ( xn , m , ) u( x, ), bmn ( x) . () () uˆ ( xn , m , ) u( x, ), bmn ( x ) The spacedomain total wavefield can be written as u ( x, ) u ( ) ( x, m , ) u ( ) ( x, m , ) . m
(B8.4)
(B8.5)
(B8.6)
m
Using the relation between the LEF and the LCB/LSB beamlets in (B8.2), the LEF () () decomposition coefficients uˆmn and uˆmn in equation (B5) 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
(B8.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 REVERSETIME CALCULATION AND LOCAL ONEWAY PROPAGATION 9.1. Background The fullwave reversetime migration method does not have the angle limitation of oneway propagators and can be used to image structures in complex areas especially for steeply dipping structures where turning waves are involved. Although the reversetime 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 reversetime migration is that the wideangle capability of fullwave propagator can generate spurious cross correlations from diving waves and backscattered 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, highpass 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 oneway 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 orderofmagnitude 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 fullwave 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 oneway 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 oneway propagators, we can extrapolate wavefield in forward or backword directions.
Figure 9.1. A simple twolayer velocity model. The horizontal and vertical lines indicate the locations where wavefields are needed for oneway extrapolations. Taking the horizontal line as an example, the wavefield u ξ and u ξ z are obtained from the fullwave calculation. Substituting these values into equation (9.1) and using oneway Green function G z r, ξ and G z r, ξ z , where superscript z denotes the oneway 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 oneway propagator only deals with narrowangle propagation and the wavefield is extrapolated only for a very short distance. The accuracy of the oneway propagator will not seriously affect the result.
101
Figure 9.2. Wavefield snapshots calculated using fullwave finitedifference 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 twolayer model. Figure 9.1 shows the velocity model. The size of the model is 451150 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 reversetime migration are calculated using a fourthorder fullwave finitedifference code. To focus our attention to the effect of backscattered waves, We mute the first arrival from the shot record. Shown in Figure 9.2 is migrated wavefield using the fullwave 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 oneway 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 backscattered 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 reversetime images for the twolayer model, where 4a is calculated using the conventional zerolag 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 oneway 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 finitedifference 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 twolayer model and major phases from different part of the structure are indicated with rays. Figure 9.6 shows the wavefield from finitedifference 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, wideangle 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 zerolag 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 oneway 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 oneway 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 finitedifference result shown in Figure 9.2, the undesired reflections have been eliminated.
104
Figure 9.4. Reversetime images for the twolayer model, (a) using conventional zerolag image condition, (b) using the image condition given by Yoon et al. (2004), and (c) using fullwave propagator plus local oneway propagator.
Figure 9.5. A simple salt dome model with steep flanks. Overlapped is the snapshot generated by finitedifference 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 fullwave reversetime calculation and local oneway propagation. Both the source and receiver waves are extrapolated using the fullwave finitedifference propagator. There is no angle limitation for these wavefields. Then the local oneway propagator is used to separate the waves into different directions. Finally, the image is calculated using directional wavefields. Using the oneway 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 orderofmagnitude smaller than if the entire wavefield is dumped out). The additional computation for the oneway 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 finitedifference 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 oneway 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 multidirectional wavefileds.
106
Figure 9.7. Depth migration image for a simple sand dome model, (a) using the conventional zerolag imaging condition, (b) using the image condition of Yoon et al. (2004), (c) using us z r u g z r u x r ux 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 diffractiontomography 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 planewave components is independent of each other, and the deconvolution filter (deconfilter) 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 crosscoupled 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 oneway 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 spacedependent deconfilter 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 AngleDomain and the Local Image Matrix The standard imaging condition in the space domain is in the form of crosscorrelation, 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 singlefrequency wavefield. In migration imaging, the final image is obtained by summing up images of all the frequencies, a doublefocusing in timedomain. In diffraction tomography, both singlefrequency tomography and multifrequency 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 anglespectra of an image field, imaging condition has been extended from space domain to the spaceangledomain (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 spaceangle 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 spacedomain is intractable. We will formulate the process of deconfiltering 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 wavenumberdomain 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 incidentscattering 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 multifrequency 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 nonuniform 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 wavenumberdomain 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 highfrequency 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 deconfilter is a division of (10.13) by (10.15): O(x, K) = I (x, k ^ )/ Â(k^ ; x) (10.17) The local image in the spacedomain 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 deconfiltering to the final multifrequency 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 singlefrequency object spectra. For a singlefrequency, 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 limitedaperture 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 deconfilter in (10.15) for a singlefrequency 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 spacedomain, 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 angledomain.
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 3D modeling project: 2nd update: The Leading Edge, 13, 949–952. Aminzadeh, F., N. Burkhard, T. Kunz, L. Nicoletis, and F. Rocca, 1995, 3D 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, Threedimensional 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. SvayLucas, 2002, Insights into migration in the angle domain: 72nd Annual International Meeting, SEG, Expanded Abstracts, 11881192. Auscher, P., 1994, Remarks on the local Fourier bases, in J.J. Benedetto and M.W. Frazier, eds., Wavelets, Mathematics and applications: CRC Press, 203218. Backus, G. and F. Gilbert, 1968, The resolving power of gross Earth data, Geophys. J. Royal Astr. Soc. 13, 247276. 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, 937954. 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, 911917. Beylkin G., Oristaglio M. and Miller D. 1985, Spatial resolution of migration algorithms, in Proceedings of the 14th International Symposium on Acoustic Imaging, pp. 155167. 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, 99108. Beylkin, G., Otistaglio, M. and Miller, D., 1986, Spatial resolution of migration algorithms, Acoustic imaging, v. 14, 155168. Billette, F. J., and S. BrandsbergDahl, 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, 12841287. Biondi, B., 2007, Angledomain commonimage gathers from anisotropic migration: Geophysics, 72, S81S91.
117
Biondi, B., and W. Symes, 2004, Angledomain commonimage gathers for migration velocity analysis by wavefieldcontinuation imaging: Geophysics, 69, 1283–1298. Bleistein, N., J. K. Cohen, and F. G. Hagin, 1987, Two and onehalf 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 Superresolution in timereversal acoustics, J. Acoust. Soc. Am. 111, 230248. BrandsbergDahl, 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, 19461949. 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, 24992503. Cao, J., and R.S. Wu, 2008, Amplitude compensation for oneway wave propagators in inhomogeneous media and its application to seismic imaging: Communications in computational physics, 3, 203221. 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, Localangle domain illumination for fullwave propagators, 78th Annual International Meeting, SEG, Expanded Abstracts, submitted. Chen, J. and Gerard T. Schuster, G., 1999, Resolution limits of migrated images, Geophysics 64, 10461053. Chen, J. and Gerard T. Schuster, G., 2001, On "Resolution limits of migrated images" (J. Chen and G.T. Schuster, Geophysics, 64, 10461053). Reply Geophys. 66, 693694 . Chen, L., R.S. Wu, and Y. Chen, 2006, Targetoriented 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 BornWKBJ inversion method for acoustic reflection data: Geophysics, 46, 15591567. 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, 259261. Collins, M.D., 1989, Applications and timedomain solution of higherorder parabolic equations in underwater acoustics, Journal of the Acoustical Society of America, 86, 10971102.
118
Collins, M. D., and E. K. Westwood, 1991, A higherorder energyconserving parabolic equation for rangedependent ocean depth, sound speed, and density: Journal of the Acoustical Society of America, 89, 10681175. Daubechies, I., S. Jaffard, and J.L. Journé, 1991, A simple Wilson orthonormal basis with exponential decay: SIAM Journal of Mathematical Analysis, 22, 554573. de Bruin, C.G.M., C.P.A. Wapenaar, and A.J. Berkhout, 1990, Angledependent reflectivity by means of prestack migration: Geophysics, 55, 12231234. De Hoop, M., J. Le Rousseau, and R. S. Wu, 2000, Generalization of the phasescreen approximation for the scattering of acoustic waves: Wave Motion, 31, 4370. Devaney, A. J., 1981, Inversescattering theory within the Rytov approximation: Optics Letters, 6, 374376. Devaney, A. J., 1982, A filtered back propagation algorithm for diffraction tomography: Ultrasonic Imag., 4, 336350. Devaney, A. J., 1984, Geophysical diffraction tomography: Trans., Inst. Electr. Electron. Eng., GE22, 313. 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, 18701873. 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, 20492052. Gelius, L.J., I., Lecomte, and H., Tabti, 2002, Analysis of the resolution function in seismic prestack depth imaging: Geophysical Prospecting, 50, 505515. Gibson, R.L. and Tzimeas, C., 2002, Quantitative measures of image resolution for seismic survey design, Geophysics, 67, 18441852. Grimbergen, J. L. T., F. J. Dessing, and K. Wapenaar, 1998, Modal expansion of one–way operators in laterally varying media: Geophysics, 63, 9951005. Han, Q., and R. S. Wu, 2005, Oneway dualdomain propagators for scalar qPwave in VTI media: Geophysics, 70, D9D17. 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. GE25, No. 4, 448455. He, Y.F., and R.S. Wu, 2007, Onereturn 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 PSwaves 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, 18181821.
119
Huang, L. J., M. Fehler, and R. S. Wu, 1999a, Extended local Born Fourier migration method: Geophysics, 64, 15241534. Huang, L. J., M. C. Fehler, P. M. Roberts, and C.C. Burch, 1999b, Extended local Rytov Fourier migration method: Geophysics, 64, 15351545. Huang, L., and M. C. Fehler, 2000, Globally optimized Fourier finitedifference migration method: 70th Annual International Meeting, SEG, Expanded Abstracts, 802805. Hubral, P., M. Tygel, and H. Zien, 1991, Threedimensional, trueamplitude zerooffset migration: Geophysics, 56, 1826. Jia, X.F., and R.S. Wu, 2007, Imaging steep salt flanks by superwide angle oneway method: 77th Annual International Meeting, SEG, Expanded Abstracts, 22652269. Jin, S., R. S. Wu, and C. Peng, 1998, Prestack depth migration using a hybrid pseudoscreen propagator: 68th Annual International Meeting, SEG, Expanded Abstracts, 18191822. Jin, S., C. C. Mosher, and R. S. Wu, 2002, Offsetdomain pseudoscreen prestack depth migration: Geophysics, 67, 18951902. Jin, S., and D., Walraven, 2003, Wave equation GSP prestack depth migration and illumination: The Leading Edge, 22, 606610. Jin, S., M., Luo, S., Xu, and D., Walraven, 2006, Illumination amplitude correction with beamlet migration, The Leading Edge, 25, 10461050. Jin, S., S. Xu, and D. Walraven, 2006, Onereturn 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 9point, finitedifference, frequencyspace, 2D scalar wave extrapolator: Geophysics, 61, 529–537. Kiyashchenko, D., R.E. Plessix, B. Kashtan and V. Troyan, 2005, Improved amplitude multioneway modeling method: Wave Motion, 43, 99115. Lambar´e, G., St´ephane Opertoz, Pascal Podvin, and Philippe Thierry, 2003, 3D ray+Born migration/inversion—Part 1: Theory, Geophysics, 68, 13481356. 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, 15381550. 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 timelapse properties of reservoirs, 73rd Ann. Internat. Mtg.: Soc. of Expl. Geophys., 15251528. Levin, S.A., 1998, Resolution in seismic imaging: Is it all a matter of perspective? Geophys., 63, 743749. Lee, M., and S.Y. Suh, 1985, Optimization of oneway 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, 985988. Luo, M.Q, J. Cao, X.B. Xie, and R.S. Wu, 2004, Comparison of illumination analyses using oneway and fullwave 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, 945948. Luo, M.Q., R.S. Wu, and X.B. Xie, 2005, True amplitude oneway propagators implemented with localized corrections on beamlets: 75th Annual International Meeting, SEG, Expanded Abstracts, 19661969. Lysmer, J., and L. A. Drake, 1972, A finiteelement 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, 22352239. 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 finitedifference and finiteelement 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 frequencydomain elastic wave modeling using weightedaveraging difference operators: Geophysics, 65, 884895. Min, D.J. and C. Shin, 2006, Refraction tomography using a waveforminversion backpropagation technique, Geophysics, 71, R21R30.. Morse, P.M., and H. Feshbach, 1953, Methods of theoretical physics: McGrawHill 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, 13791382 Muerdter, D., and D. Ratcliff, 2001a, Understanding subsalt illumination through raytrace modeling, Part 1: Simple 2D salt models: The Leading Edge, 20, 578594.
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, 688687. Muerdter, D., and D. Ratcliff, 2001c, Understanding subsalt illumination through raytrace modeling, Part 3: Salt ridges and furrows, and impact of acquisition orientation: The Leading Edge, 20, 803816. Mufti, I.R. J.A. Pita, and R.W. Huntley, 1996, Finitedifference depth migration of explorationscale 3D seismic data: Geophysics, 61, 776794. Mulder, W.A., and R.E. Plessix, 2004, A comparison between oneway and twoway waveequation migration: Geophysics, 69, 14911504. 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, 13571370. Papanicolaou, G., L. Ryzhik, and K Solna, Statistical stability in timereversal, SIAM J. Appl. Math., 64, 11331155. Pratt, R.G., 1999, Seismic waveform inversion in the frequency domain, Part 1: Thoery and verification in a physical scale model: Geophysics, 64, 888901. Pratt, R.G. and Worthington, M.H., 1988, The application of diffraction tomography to crosshole seismic data: Geophysics, 53, 12841294. ———, 1990, Inverse theory applied to mulitplesource crosshole tomography I: Acoustic waveequation method: Geophys. Prosp., 38, 287310. Priolo, E. and Chiaruttini, C., 2003, Analytical and numerical analysis of tomographic resolution with bandlimited signals, Geophysics, 68, 600613. Prucha, B. B., and W. Symes, 1999, Angledomain commonimage gathers by waveequation migration: 69th Annual International Meeting, SEG, Expanded Abstracts, 824–827. Qian, S. and D. Chen, 1996, Joint TimeFrequency Analysis, Methods and Applications: Prentice Hall PTR. Rickett, J.E., 2003, Illuminationbased normalization for waveequation depth migration, geophysics, 68, 13711379. Rickett, J., and P. Sava, 2001, Offset and angle domain commonimage gathers for shotprofile migration: 71st Annual International Meeting, SEG, Expanded Abstracts, 1115–1118. Rickett, J.E. and P.C. Sava, 2002, Offset and angledomain common imagepoint gathers for shotprofile migration: Geophysics, 67, 883889. Ristow, D., and T. Rühl, 1994, Fourier finitedifference migration: Geophysics, 59, 1882–1893. Rosales, D.A., S. Fomel, B. Biondi, and P. Sava, 2008, Waveequation angledomain commonimage gathers for converted waves: Geophysics, 73, S17S26. Sava, P. C., and S. Fomel, 2003, Angledomain commonimage gathers by wavefield continuation methods: Geophysics, 68, 1065–1074.
122
Schneider,W. A., and G. A.Winbow, 1999, Efficient and accurate modeling of 3D seismic illumination: 69th Annual International Meeting, SEG, Expanded Abstracts, 633–636. Schultz, P. S., and J. F. Claerbout, 1978, Velocity estimation and downwardcontinuation by wavefront synthesis: Geophysics, 43, 691–714. Schuster, G.T., and Hu, J., 2000, Green’s function for migration: Continuous recording geometry, Geophysics, 65, 167175. Sheng, J., A. Leeds, M. Buddensiek, and G.T. Schuster, 2006, Early arrival waveform tomography on nearsurface refraction data, Geophysics, 71, U47U57. Shin, C., and H. Sohn, 1998, A frequencyspace 2D scalar wave extrapolator using extended 25point finitedifference 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., 10551058. Sjoeberg, T.A. and Lecomte, I., 2003, 2D deconvolution of seismic image blur, 73rd Annual International Meeting, SEG, Expanded Abstracts, 10551058. Snieder, R. and Lomax, A., 1996, Wavefield smoothing and the effect of rough velocity perturbations on arrival times and amplitudes, Geophys. J. Int. 125, 796812. Spetzler, J. and Snieder, R., 2001, The effects of smallscale heterogeneity on the arrival time of waves, Geophys. J. Int. 145, 786796. 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 finitedifference method, Geophys. J. Int. 148, 165178. Spetzler, J. and Snieder, R., 2004, Tutorial: The Fresnel volume and transmitted waves, Geophysics, 69, 653663. Stekl, I., and R. G. Pratt, 1998, Accurate viscoelastic modeling by frequencydomain finite differences using rotated operators: Geophysics, 63, 1779–1794. Steinberg, B. Z., 1993, Evolution of local spectra in smoothly varying nonhomogeneous environmentslocal 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 reversetime depth migration of prestack elastic seismic data: Geophysics, 66, 15191527. Tabti, H., L. Gelius, and I. Lecomte, 2001, Discussion On “Resolution limits of migrated images” (J. Chen and G. T. Schuster, Geophysics, 64, 10461053). 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 wideangle wavefield extrapolators and scattering operators: Geophysical Journal International, 163, 308323. 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, 918931. Wang, Y. and R.S. Wu, 2002, Beamlet prestack depth migration using local cosine basis propagator: 72nd Annual International Meeting, SEG, Expanded Abstracts, 13401343. Wang, Y., R. Cook, and R.S. Wu, 2003, 3D local cosine beamlet propagator: 73rd Annual International Meeting, SEG, Expanded Abstracts, 981984. 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, 22612264. Wapenaar, C. P. A., and J. L. T. Grimbergen, 1996, Reciprocity theorems for oneway 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, 202207. Woodward, M., 1992, Waveequation tomography: Geophysics, 57, 1526. Wu, R. S., 1994, Wideangle elastic wave oneway propagator in heterogeneous media and an elastic wave complexscreen method: Journal of Geophysical Research, 99, 751766. Wu, R. S., 1996, Synthetic seismograms in heterogeneous media by onereturn approximation: Pure and Applied Geophysics, 148, 155173. Wu, R.S., F.V. Araujo and L.J. Huang, 1994, Multifrequency backscattering tomography for constant and vertically varying backgrounds, Int. J. of Imaging Syst. & Tech., 5, 721. 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 GaborDaubechies frame propagator: 63rd Annual International Meeting, EAGE, Expanded Abstracts, 7478. Wu, R.S., and L. Chen, 2002, Mapping directional illumination and acquisitionaperture efficacy by beamlet propagators: 72nd Annual International Meeting, SEG, Expanded Abstracts, 13521355. Wu, R.S. and Chen, L., 2002b, Wave propagation and imaging using GaborDaubechies beamlets: "Theoretical and Computational Acoustics", World Scientific, New Jersey, 661670.
124
Wu, R.S., and L. Chen, 2003, Prestack depth migration in angledomain 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, S147S159. Wu, R.S., and L. Chen, 2006b, Targetoriented beamlet migration based on GaborDaubechies 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 dipresponse, 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 superwide angle oneway waves by wavefront reconstruction: 76th Annual International Meeting, SEG, Expanded Abstracts, 29762980. Wu, R.S., M.Q. Luo, S.C. Chen, and X.B. Xie, 2004, Acquisition aperture correction in angledomain and trueamplitude imaging for wave equation migration: 74th Annual International Meeting, SEG, Expanded Abstracts, 937940. Wu, R.S., and J. Mao, 2007, Beamlet migration using local harmonic bases: 77th Annual International Meeting, SEG, Expanded Abstracts, 22302234. Wu, R.S. and Toksöz, M.N., 1987, Diffraction tomography and multisource holography applied to seismic imaging, Geophysics, 52, 1125. 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, 10081011. Wu, R.S., Y. Wang and M.Q. Luo, 2003, Localcosine beamlet migration for 3D complex structures: 8th International Congress of the Brazilian Geophysical Society, 1418. 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 higherorder, finitedifference schemes in 3D reversetime migration: Geophysics, 61, 845856. Wu, R.S. ,Wu, Y.L., Liu, C.R., Peng, C.Y., and M.J. Wang, 1977, Multifrequency synthetic detecting holography, Acta Geophysica Sinica, 20, 3338. Xie, X. B., Z., Ge and T., Lay, 2005, Investigating explosion source energy partitioning and Lgwave excitation using a finiteDifference plus slowness analysis method, Bull. Seism. Soc. Am. 95, 24122427.
125
Xie, X. B., S.W. Jin, and R. S. Wu, 2003, Threedimensional 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, 933936. Xie, X.B., S.W. Jin and R.S. Wu, 2006, Waveequationbased seismic illumination analysis: Geophysics, 71, S169S177. 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, 18111814. Xie, X. B., and R. S. Wu, 2005, Multicomponent prestack depth migration using elastic screen method: Geophysics, 70, S30S37. Xie, X.B., and R.S. Wu, 2002, Extracting angle domain information from migrated wave field: 72nd Annual International Meeting, SEG, Expanded Abstracts, 13601363. Xie, X.B., R.S. Wu, M. Fehler and L. Huang, 2005, HSeismic resolution and illumination: A waveequationbased analysisH, 75th Annual International Meeting, SEG, Expanded Abstracts, 18621865. Xie, X.B., and R.S., Wu, 2006, A depth migration method based on the fullwave reversetime calculation and local oneway propagation, 76th Annual International Meeting, SEG, Expanded Abstracts, 23332337. Xie, X.B., Wu, R.S., M.Fehler and L. Huang, 2005, Seismic resolution and illumination: A waveequationbased analysis: 75th Annual International Meeting, SEG, Expanded abstracts, 18621865. Xie, X.B., and H., Yang, 2008, A fullwave 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, Commonangle migration: A strategy for imaging complex media: Geophysics, 66, 18771894. 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, 3D prestack fullwavefield inversion: Geophysics, 60, 18051818. Yang, H., X.B. Xie, M.Q. Luo, and S.W. Jin, 2008, Target oriented fullwave 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 reversetime migration using the acoustic wave equation: An experience with the SEG/EAGE data set: The Leading Edge, 22, 3841.
126
Yoon, K., K.J. Marfurt, and W. Starr, 2004, Challenges in reversetime migration: 74th Annual International Meeting, SEG, Expanded Abstracts, 10571060. Yu, J., and Schuster, G.T., 2003, 3D prestack migration decomvolution: 73rd Annual International Meeting, SEG, Expanded Abstracts, 16511654. Zhou, C., W. Cai, Y. Luo, G.T. Schuster, and S. Hassanzadeh, 1995, Acoustic waveequation traveltime and waveform inversion of crosshole seismic data: Geophysics, 60, 765773. Zhou, C., G.T. Schuster, S. Hassanzadeh, and J.M. Harris, 1997, Elastic waveequation traveltime and waveform inversion of crosswell seismic data: Geophysics, 62, 853868. 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 oneway wave equations: Inverse Problems, 19, 11131138. Zhang, Y., G. Zhang, and N. Bleistein, 2005, Theory of trueamplitude oneway wave equations and trueamplitude commonshot migration: Geophysics, 70, E1–E10. Zhang, Y., S. Xu, and G.Q. Zhang, 2006, Imaging complex salt bodies with turningwave oneway wave equation: 76th Annual International Meeting, SEG, Expanded Abstracts, 2323–2326.
127