Seismic Illumination Analysis with One-Way Wave Propagators Coupled with Reflection/Transmission Coefficients in 3D Complex Media

Seismic illumination analysis can help us to better understand how various seismic acquisition parameters and configuration affect seismic image qualities. It allows us to design effective seismic survey systems to provide reliable high-resolution seismic images for complex structures in both theoretical and exploration seismology. It is also a useful tool for estimation of the potential detecting power of a specific acquisition system for a given velocity structure of the medium. Most existing techniques for calculating illumination intensity distributions can be divided into two categories: one is based on ray tracing methods (Berkhout, 1997; Schneider, 1999; Bear et al., 2000; Muerdter et al., 2001a, 2001b, 20001c); and the other based on wave equation methods (Wu & Chen, 2002; Xie et al., 2006; Wu & Chen, 2006; Sun et al., 2007). The ray tracing is quite efficient and inexpensive. However, it bears large errors in complex areas because of the multi-path problem, the high frequency approximation and the singularity problem of the ray theory (Hoffmann, 2001). These factors limit its accuracy severely in complex areas. The wave-equation-based methods can provide reliable and accurate illumination intensity distributions in complex media. Full-wave finite-difference (FD) method is commonly employed for seismic forward modelling, i.e., wave propagation simulation for a known interval velocity model. But it is too expensive for illumination analysis in industrial application, especially for 3D case. Recently, One-way wave propagators have been widely used to seismic modelling and migration for its huge memory-saving and high computational efficiency. Localized illumination methods based on the Gabor-Daubechies frame decomposition or the local slant stack have been presented, which can give local information both in space and direction (e.g., Xie and Wu, 2002; Wu et al., 2000; Xie et al., 2003; Wu & Chen, 2006; Mao & Wu, 2007; Cao & Wu, 2008). However, its computational burden is still unacceptable due to the wavefield decomposition in every spatial point.


Introduction
Seismic illumination analysis can help us to better understand how various seismic acquisition parameters and configuration affect seismic image qualities. It allows us to design effective seismic survey systems to provide reliable high-resolution seismic images for complex structures in both theoretical and exploration seismology. It is also a useful tool for estimation of the potential detecting power of a specific acquisition system for a given velocity structure of the medium. Most existing techniques for calculating illumination intensity distributions can be divided into two categories: one is based on ray tracing methods (Berkhout, 1997;Schneider, 1999;Bear et al., 2000;Muerdter et al., 2001aMuerdter et al., , 2001b; and the other based on wave equation methods (Wu & Chen, 2002;Xie et al., 2006;Wu & Chen, 2006;Sun et al., 2007). The ray tracing is quite efficient and inexpensive. However, it bears large errors in complex areas because of the multi-path problem, the high frequency approximation and the singularity problem of the ray theory (Hoffmann, 2001). These factors limit its accuracy severely in complex areas. The wave-equation-based methods can provide reliable and accurate illumination intensity distributions in complex media. Full-wave finite-difference (FD) method is commonly employed for seismic forward modelling, i.e., wave propagation simulation for a known interval velocity model. But it is too expensive for illumination analysis in industrial application, especially for 3D case. Recently, One-way wave propagators have been widely used to seismic modelling and migration for its huge memory-saving and high computational efficiency. Localized illumination methods based on the Gabor-Daubechies frame decomposition or the local slant stack have been presented, which can give local information both in space and direction (e.g., Xie and Wu, 2002;Wu et al., 2000;Xie et al., 2003;Wu & Chen, 2006;Mao & Wu, 2007;Cao & Wu, 2008). However, its computational burden is still unacceptable due to the wavefield decomposition in every spatial point.
In this chapter, we developed a 3D one-way wave-equation-based illumination method for complex media. It combines the one-way propagators coupled with reflection/transmission (R/T) coefficients and the phase encoding techniques. The one-way propagators coupled with R/T coefficients were firstly developed by Sun et al. (2009) for the 2D case. The R/T operators not only account for amplitude variations with incident angles across interfaces, but also accommodate to complex media with steep dip angle and large lateral velocity contrast. Firstly, we extended the method to the 3D case in complex structures starting from generalized Lippmann-Schwinger equations and applied the method to calculate illumination intensity distributions. Then, the phase encoding technique in frequency domain (Romero et al., 2000) was employed to reduce the computational cost further by calculating a number of shots together, which can apply to any frequency domain illumination methods. In the following sections, the accuracy of the R/T propagators, defined by the corresponding dispersion relationship, was analyzed by comparison with the exact solution and other popular one-way propagators known as the split-step Fourier (SSF) method (Stoffa et al., 1990) and the generalized screen propagator (GSP) (de Hoop et al., 2000). Then, special attention was given to the implementation of numerical procedures to improve efficient computation efficiency and huge memory savings. Several numerical examples were given to show its capabilities to handle complex media. As the actual interval velocity model or geological model cannot be known exactly, the effects of velocity error on illumination intensity distributions were evaluated by numerical examples. Finally, we showed some practical applications of target-oriented illumination analysis in seismic acquisition and survey design, which is fundamental to exploration geophysics.

Methodology
I n t h i s s e c t i o n , w e f i r s t d e r i v e 3 D v e rsion of one-way propagators coupled with reflection/transmission (R/T) coefficients. Then, we briefly introduce the definition of seismic illumination and theories of phase encoding. The accuracy of R/T propagators is analyzed by comparison with the exact solution and other popular one-way propagators. Finally, we present some applications of seismic illumination.

One-way Lippmann-Schwinger wave equation
Generally, most media can be sliced into heterogeneous slabs perpendicular to the major propagation direction. Figure 1 depicts the geometry of such a heterogeneous slab denoted by  with the top interface 0  , the bottom interface 1  , and the thickness z  . The velocity distribution in the slab is denoted by () v r where r is the position vector, and its reference velocity is 0 v . We start with the scalar Helmholtz equation for a time-harmonic wavefield where / n denotes differentiation with respect to the outward normal of the boundary  .   (3) and (4) into (2) and considering the "boundary naturalization" of the integral equations, that is, a limit analysis when the "observation point" r approaches the boundary  and tends to coincide with the "scattering point"   r , we obtain the following generalized Lipmann-Schwinger integral equation www.intechopen.com for all     r , where the coefficient () 1 / 2 C  r for a flat  . This is a wave integral equation that is equivalent to the Helmholtz equation (1) and describes two-way wave propagation in the heterogeneous slab. It is a Fredholm integral equation of the second kind with the existence and uniqueness of its solution for both interior and exterior Helmholtz problems assured by the classical Fredholm's theorems of integral equations. Numerical studies of this equation based on the discrete-wavenumber boundary methods have been conducted for wave propagation simulation (Fu and Bouchon, 2004) in piecewise heterogeneous media that the earth presents. Because numerous matrix operations are involved and the matrix for each frequency -component computation must be inverted, the numerical methods are computationally intensive at high frequencies.
Equation (5) describes wave propagation in the space-frequency domain. Formulating it in the frequency-wavenumber domain using a plane-wave expansion will lead to quite different numerical schemes. In Figure 1, we assume wave propagation along the z-axis, crossing the slab from the slab entrance 0  to the exit 1  . Let ( where ( , ) This is actually the Born approximation applied to the slab. It implies that the heterogeneity of the slab is represented by its top/bottom interfaces and consequently requires the slab is thin enough with respect to the wavelength of incident waves. Substituting equations (6), (7), and (8) into (5) and noting that each inner integral is a Fourier transforms, we obtain Equation (9) is a wavenumber-domain wave equation that describes two-way wave propagation in the heterogeneous slab, including multiple forward and back scatterings between 0  and 1  .
For one-way wave propagation using the matching solution techniques, further simplification should be made to Equation (5) by reducing it to one-way version. From Equation (9), we see that two-way wave propagation involves two terms: the displacement () u r and acoustic pressure gradient () q r . In practical, we do not often measure both () u r and () q r at a given level. The pressure gradient () q r at the slab entrance 0  can be dropped www.intechopen.com by choosing 0  as an acoustically soft boundary (Dirichlet boundary condition) which leading to one-return approximation shown in Figure 2. This Rayleigh-type integral representation is valid if we neglect back scatterings. With this choice, no energy returns from the upper boundary 0  and multiple reflections between 0  and 1  can be avoided.
This choice updates Equation (9) to To account for the effect of transmission and refraction at 1  in a natural manner, we need to build a boundary integral equation in the medium immediately below the slab with the radiation conditions applied to the bottom boundary of the medium, where z k is the wavenumber related to the medium immediately below 1  , and (, ) Equation (13) is one-way downward wave equation with one-return approximation that accounts for both the accumulated effect of forward scattering by volume heterogeneities inside the slab and the transmission between adjoining slabs on wave amplitude and phase. Next, we derive the formula of one-way upward waves. For upward wave propagation, it has up . We rewrite the volume integration over the slab in Equation (5) using the rectangular rule Substituting equations (6), (7) and (14) into (5), and noting that each inner integral is a Fourier transform, we have To account for the effects of the reflection and transmission at 0  and 1  on backward wave propagation, we build two boundary integral equations: Substituting the Green function to equations (16) and (17), and using the approximation (, ) (, , substituting into equation (15), we obtain This is one-way upward wave equation with one-return approximation. Equations (13) and (18) account for both the effects of forward and backward scattering by volume heterogeneities inside a slab and the R/T between adjoining slabs on wave amplitude and phase.

One-way propagators 2.2.1 First-order separation-of-variables screen propagators
We consider the following splitting operator decomposition (1 )   T k in equation (21) can be approximated by the following rational expansion where the coefficients j a and j b are independent of n . Submitting this equation into Equation (21) Coefficients in Equation (23) can be determined numerically by an optimization procedure using the least-squares method. Because of the mathematical properties and approximation behavior of rational functions (Trefethen and Halpern, 1986), equation (23) should be wellposed especially for the lower-order terms. In practice, its first-order equation or at most the second-order equation is adequate for common one-way propagation in large-contrast media with wide propagation angles in seismology.
In what follows, we formulate these separation-of-variables screen propagators (SVSP) by a Fourier-transform-based representation for numerical implementation. Substituting Equations (23) The computational time with Equation (25) is almost the same as traditional split-step Fourier solutions, but with high accuracy close to the generalized screen propagator (GSP). Note that there is a singularity in equation (25) when z k and z k approach zero simultaneously which leads to an instability of the algorithm. This can be avoided using the following relation (Huang et al., 1999): where  is a small real number. The R/T progators can be also applied to other known one-way propagators, e.g., SSF, FFD, due to the same or similar operator structure. Here, we give the common form in the frequency-wavenumber domain as The present method is a separation-of-variables Fourier marching algorithm. The whole medium is sliced into a stack of thin slabs perpendicular to the main propagation direction. The implementation procedures may be summarized as follows: 1. Slice the medium into a stack of slabs perpendicular to the propagation direction (i.e., along the z-axis direction). 2. Choose a reference velocity in each slab to make it a perturbation slab represented by the acoustic refractive index. 3. Interact with the split-step terms and Fourier transform the wave fields at the entrance of each slab into the frequency-wavenumber domain.
4. Conduct the linear interpolation and propagate the wave fields through the slab in the frequency-wavenumber domain by multiplying a phase shift using the reference velocity.

Interact with the R/T coefficients to get the transmitted fields and reflected fields and
inverse Fourier transform the wave fields at the exit of the slab into the frequency-space domain. 6. Repeat steps 3 to 5 slab by slab until the last slab.

Analysis of relative phase errors
The corresponding dispersion relation is shown in Equation (23) Figure 3a compares the angular spectra of the first-order separation-of-variables screen propagator (SVSP1) with the GSP and SSF propagators under a 5% relative phase error. We see that the SVSP1 works perfectly for all the n values and all propagation angles, almost approaching that of the GSP. To investigate the global properties of equation (30), we evaluate the first-and second-order separation-of-variables approximations across three different n values by comparing their dispersion circles with the exact, GSP, and SSF dispersion relations ( Figure 3b). As expected, the approximation can be an exact expression in the small-angle pie slice (0 ) x k  for all the n values. This accuracy comparison of the first-and second-order screen propagators with other Fourier propagators demonstrates a quick converge of the regional separation-of-variables approximation in the low-order terms.

Definition of seismic illumination
Seismic illumination stands for the distribution of seismic wave energy underground. Since seismic wave energy is located at a frequency band cantered at peak frequency, the illumination map can be obtained by summing up seismic wavefields over a small frequency band. We define multi-frequency point source illumination (,) s Iz x as 0 0 1/2 2 (,) (,; ) where  is angular frequency, 0  is peak frequency, and   is half of the frequency band.
where (,; ) Sz  x and (,; ) Gz  x are synthetic plane-wave wavefields of source and geophone wavefields at location (,) z x for a frequency, respectively.

Computational efficiency
The one-way propagators are shuttled between the space domain and the wavenumber domain via fast Fourier transform (FFT). Thus, the computational efficiency of our method is dependent on the performance of FFT. It is necessary to pick an efficient FFT for fast implementation of our method. As known, lots of FFT codes are contributed to us, which have different computation cost. Here, we test two of the fastest FFT packages, i.e., FFT in west (FFTW) (Frigo and Johnson, 2005) and Intel Math Kernel Library (MKL). Both are free packages under the terms of the GNU General Public License. Although computer technologies have been advanced greatly, people do not satisfy with computational cost yet. Thus, it is still necessary to choose FFT as fast as possible. The benchmark is performed on a personal computer. The configuration and compiler environment is shown in Table 1.
To test the performance of the two FFT packages, we take a 2D Gaussian function as input data, defined as where xi x  , yj y   , 1,2, , We compute the floating point operations per second (FLOPS). The formula is given by where NN xN y is the length of FFT, t is the CPU time of performing forward and inverse FFT, its unit is second. Noting that the FLOPS is not true one, which is a normalized value by 2 5l o g () NN , which is the redundancy of the radix-2 Cooley-Tukey algorithm.  Figure 4 shows the computational speed of the two FFT packages. The FLOPS of Intel MKL is much higher than the one of FFTW, especially at location of commonly used 2D FFT size ( 512x1024-512x8192 ). Thus, we choose Intel MKL to perform FFT in our method. www.intechopen.com

Applications
In this section, we will demonstrate the numerical and practical applications of seismic illumination in exploration geophysics. Firstly, we verify its accuracy by comparing illumination map with prestack depth migration for a 2D salt model. Secondly, we take 2D SEG/EAGE salt model as an example to obtain a criterion for seismic illumination and migration. Finally, we take a 3D fault block model in eastern China to show its powerful applications in industry. Figure 5 displays a theorectical model with a salt body at the center. The horizontal distance is 12000m and the depth is 1440m. The space sampling is 25m and 6m in the horizontal and depth directions, respectively. There are 120 receivers per shot. The shot interval is 50m and the receiver interval is 25m. The plane-wave source illumination map and prestack depth migration (PSDM) profile are shown in Figure 6. The result of illumination map and migration profile are quite similar. While it shows larger energy in the illumination map, it shows good phases and amplitudes in PSDM profile, and vice versa. In Figure 6a, there is an obvious illumination shade due to the exsitance of salt body. Accordingly, it has poorly imaged at the same position in figure 6b. Geologic structures have great influence on waves propatation udnerground. This leads to uneven distribution of seismic energy. To demonstrate the influence of geologic structures on illumination energy, we give two point source illumination maps at different location (4000m and 8000m) in Figure 7. In Figure 7a, the shot is located at leftside of salt body, where the energy distributes uniformly. However, the distribution of seimic wave energy is seriously uneven while the shot locates at the rightside of salt body (see figure 7b) and will lead to poor images.

Seismic illumination and PSDM
In this section, we clarify two criterions of seimic illumination by comparing illumination map with PSDM profile. With the help of the two criterions, people can predict the results of seismic acquisition and migration via illumination distribution. The two criterions are 1) strong illumination energy and 2) uniform distribution of illumination. Here, we take 2D SEG/EAGE salt model as an example for simplicity. The most important reason to pick 2D SEG/EAGE salt model is that the model is a standard model used to test various migration algorigthms. Figure 8a shows the 2D SEG/EAGE salt model. The horizontal and vertical sampling numbers are 1290 and 300, respectively. The spatial sampling interval is 12.5m. We use Fourier finite-difference method (FFD) (Ristow and Rühl, 1994) to obtain the post-stack migraion result shown in Figure 8b. Figure 8c shows the illumination map calculated by plane wave source. The source function is Ricker wavelet with dominant frequency 15Hz.
The frequency band is 2Hz, from 14Hz to 16Hz. The incident angle of plane wave ranges from 50   to 50  . In Figure 8c, two shade areas subsalt can be observed obviously. This is because seismic waves cannot penetrate the salt and are reflected. Comparing the migration results and the illumination map, we find that both have good consistence with each other. That means, the poorer the illumination is, the worse the migration is, and vice versa.   Figure 8a and (c) the corresponding illumination intensity map

Effects of velocity uncertainties
Same as the prestack depth migration (PSDM), the illumination maps are also calculated based on interval velocity model in depth domain. In general, the actual interval velocity model or geological model cannot be known exactly. The quality of PSDM is dependent on the accuracy of interval velocity model. Thus, it is necessary to evaluate the effects of velocity errors on illumination intensity distributions.
In this section, we take a 3D real geologic model ( Figure 9) in Eastern China as an example to investigate the influence of velocity error on the illumination intensity. The geologic structure in this region is very complex. Small faults/blocks are quite common. The migration images from this area are relatively poor. The vertical slices of the model are shown in figure 9a and 9b. The horizontal sampling numbers are 501 and 201, respectively. The depth sampling number is 209. The spatial sampling interval is 25m in all directions. The depth of target zone is 3000-4000m. The source function is Ricker wavelet with dominant frequency 15Hz. The frequency band is 2Hz, from 14Hz to 16Hz. The strike of the geologic structure is along the x-axis. The incident angle of plane wave ranges from 50   to 50  in the x direction, and from 20   to 20  in the y direction. Here, we first give an illumination map with a point source in Figure 10. In this example, we use correct interval velocity. The energies in target zone are distributed unuiformly, which illustrates complex structures have serious effects on wave propagations.   We compare the illumination maps calculated with 90%, 100% and 110% correct velocity, shown in Figure 11-13, respectively. Although the illumination maps are obtained by using incorrect velocities, their characteristics are almost same with each other. Thus, we can conclude that the small velocity error has little influences on illumination intensity.

Target-oriented illumination analysis
As the criterions stated above, we can obtain a good quality image of subsurfaces, if the illumination intensity is strong and distribute uniformly. Or else a poor quality image will be resulted. From Figure 12, the illumination in target zone is weak and does not distribute uniformly, which will lead to poor image underground.
To obtain an image with good quality, people are always trying to satisfy the two criterions. An direct way is to enhance the illumination intensity by adding shot excitations in the field.
The key of the method is where to place the shots and receivers in order to enhance the illumination in target zones. Now, we try to locate effective shot postions according to the reciprocity theory. The illumination at the surface can be considered as the contribution to the target zone while shots are excited at the surface. We firstly place the seismic sources in the target zone. Then, we propagate the wavefields to the surface and calculate the illumination intensity at the surface. Figure 14 shows us the illumination at the surface where the source are placed in the target zone. It can be found that the illumination does not distribute uniformly, since the wavepaths are distorted while waves run in the complex media. The results of Figure 14 can help us to determinate where to add shots. After the effective locations of shots excitation are obtained, we next ensure locations of corresponding receivers. Firstly, we put the source at where has largest contributions to improve the illumination intensities of target zones. Then, the reflectivities in the complex media are set to zero except in the target zones. According to equations (13) and (18), the wavefields are propagated in the complex media. The wavefields recorded at the surface are reflected only from the target zone. Thus, we can use the target-reflected illumination intensities to locate where to receive the wavefields. Figure 14 show the optimal receivers position for a given source at (6200,1500)m. In Figure 13 and 14, both optimal sources and receivers positions are irregular, since the wavefields are distorted while running in the compelx media. Thus, target-oriented illumination analysis would totally change the conventional way of designing seismic acquisition system.

Conclusion
Starting from generalized Lippmann-Schwinger equations, we developed a 3D one-way wave-equation-based illumination method in complex media. It combines the one-way propagators coupled with reflection/transmission (R/T) coefficients and the phase encoding techniques. The R/T operators not only account for amplitude variations with incident angles across interfaces, but also accommodate to complex media with steep dip angle and large lateral velocity contrast.
Several numerical examples are given to illustrate its resolving capabilities for complex media. In this chapter, we have demonstrated the numerical and practical applications of seismic illumination in exploration geophysics. Numerical examples shows that the illumination maps and results of post-and pre-stack depth migration are consistent with each other. Two criterions i.e. strong illumination energy and uniform distribution of illumination are obtained by comparing the illumination intensity map with the prestack depth migration for 2D SEG/EAGE salt model. The velocity errors to the illumination intensities are investigated, which shows the velocity error has little influences on illuminations. The target-oriented illumination analysis has been applied to design seismic acquisition, which is fundamental to the exploration geophysics.

Acknowledgment
The helpful discussions with T. Jiang and X. J. Wan are greatly appreciated. This work was supported by China Postdoctoral Science Foundation (Project No. 20100480447).