Wave equation depth migration using complex Padé approximation

Wave equation depth migration using complex Padé approximation

Reynam C. Pestana^I; Jacira C.B. Freitas^II; Jessé C. Costa^III

^IInstitute of Physics & Research Center in Geophysics and Geology, Federal University of Bahia, Salvador, BA, Brazil. Phone: [phone]; Fax: [phone] - E-mail: [e-mail]

^IIInstitute of Physics & Research Center in Geophysics and Geology, Federal University of Bahia, Salvador, BA, Brazil. Phone: [phone] Extension no.:6614; Fax: [phone] - E-mail: [e-mail]

^IIIFaculty of Geophysics, Federal University of Pará, Guamá, 66050-970 Belém, PA, Brazil. Phone: [phone] Extension no.: 41 - E-mail: [e-mail]

ABSTRACT

We propose a new method of depth migration based on a constant density variable velocity wave equation in the space-frequency domain. A complex Padé approximation of the wave equation evolution operator is used for wavefield extrapolation. This method mitigates the inaccuracies and instabilities due to evanescent waves and produces images with fewer numerical artifacts than those obtained with a real Padé approximation of the exponential operator, mainly in media with strong velocity variations. Tests on zero-offset data from the SEG/EAGE salt model and the 2D Marmousi prestack dataset show that the proposed migration method can handle strong lateral variations and also has a good steep dip response. We compare the results of the proposed method with those obtained using split-step Fourier (SSF), phase shift plus interpolation (PSPI) and Fourier finite-difference (FFD) methods.

Keywords: depth migration, one-wave equation, complex Padé approximation, prestack migration.

RESUMO

Propomos um novo método de migração em profundidade baseado na solução da equação da onda com densidade constante no domínio da freqüência. Uma aproximação de Padé complexa é usada para aproximar o operador de evolução aplicado na extrapolação do campo de ondas. Esse método reduz as imprecisões e instabilidades devido às ondas evanescentes e produz imagens com menos ruídos numéricos que aquelas obtidas usando-se a aproximação de Padé real para o operador exponencial, principalmente em meios com fortes variações de velocidades. Testes em dados de afastamento nulo do modelo de sal SEG/EAGE e nos dados de tiro comum 2-D Marmousi foram realizados. Os resultados obtidos mostram que o método de migração proposto consegue lidar com fortes variações laterais e também tem uma boa resposta para refletores com mergulhos íngremes. Os resultados foram comparados àqueles resultados obtidos com os métodos split-step Fourier (SSF), phase shift plus interpolarion (PSPI) e Fourier diferenças-finitas (FFD).

Palavras-chave: migração em profundidade, equação unidirecional da onda, aproximação de Padé complexa, migração pré-empilhamento.

INTRODUCTION

Wave equation migration methods perform better imaging than ray-based migration methods in media with strong lateral velocity variations. The majority of wave equation methods are based on solving one-way acoustic wave equation. There are many different methods to numerically solve a given one-way wave equation. These methods can be grouped into three classes: Fourier methods (Stoffa et al., 1990; Gazdag & Sguazzero, 1984) solved in the wave number domain; finite-difference methods (FD) (Claerbout, 1985; Hale, 1991) and mixed methods that are a linear combination of spectral and Fourier finite-difference methods (FFD) (Ristow & Rühl, 1994). In both FD and FFD migration methods, the purely real Padé approximation or a simple Taylor expansion is usually used to approximate the square-root operator. Real operators cannot handle evanescent waves (Millinazzo et al., 1997), to deal with this problem they proposed a complex Padé approximation.

The complex Padé expansion was previously used in applied geophysics. Zhang et al. (2003) used the method in finite-difference migration however its implementation was not efficient for wide-angle migration. Zhang et al. (2004) derived a split-step complex Padé-Fourier solution for one-way equation using both [1/1] and [2/2] Padé approximations. Recently, Amazonas et al. (2007) used the complex Padé approximation for wide-angle FD and FFD prestack depth migration. These works show that complex Padé approximation mitigates the inaccuracy and instability due to evanescent waves.

We use the exact acoustic wave equation to derive the downward continuation operator, instead of starting from an one-way wave equation. The propagation of the wavefield from one depth level to the next is carried out by the complex Padé approximation of the exact downward continuation operator. The discretized approximation of this operator produces a tridiagonal linear system which is solved to obtain the extrapolated wavefield from the previous one.

In the following sections, we derive the FD Padé solution for the wave equation using the complex Padé approximation. Then we present zero-offset migration results for the SEG/EAGE salt model. Next, the 2D Marmousi dataset is used to test the prestack shot-record depth migration scheme that is developed, based on the FD method using complex Padé approximation propagators. The cross-correlation between source and receiver wavefields is applied for prestack imaging condition. A comparison with the standard (real-valued) FD, FFD, Fourier split-step, PSPI and our method is also given. The test results show that our FD migration method is capable of imaging structures with strong lateral velocity variation.

THEORY AND THE METHOD

The finite-difference solution method is based on a depth extrapolation of the temporally transformed constant density acoustic wave equation given by:

where P (x, y, z, ω) denotes the Fourier transform of the pressure field p (x, y, z, t), ω is the temporal frequency, x and y are the horizontal coordinates, z is the depth and c(x, y, z) is the velocity.

The wave equation may be spatially discretized on a uniform mesh (x, y, z) followed by a selection of a second derivative approximation to transform it into a system of ordinary differential equations (Kosloff et al., 2006) as follows,

where,

and c_o is a constant velocity (usually chosen as the minimum velocity in each layer).

We assume that within each layer (z, z + Δz) the velocity is vertically constant. Then the upward propagating solution to (2) within a layer can be formally written as

Introducing , we have

The operator D in (3) contains both positive and negative values which respectively correspond to propagating and evanescent waves (Kosloff et al., 2006).

A complex treatment of the exponential term is applied to mitigate inaccuracy and instabilities due to the evanescent waves (Milinazzo et al., 1997; Zhang et al., 2004). This can be done using a complex Padé approximation of the square root.

Before applying the Padé approximation, the branch cut of the square-root is rotated by an angle θ. The operator in equation (5) can be rewritten as

where δ' = δ e^{i θ/2} and D' = D e^{-i θ} - 1 . Then

where, Z = - 1 .

Applying Padé approximation to we have (Zhang, et al., 2003):

where N defines the order of the Padé expansion, a_n and b_n are real-valued coefficients. This expansion permits us to write

Considering only the first-order approximation, equation (8) becomes

where a = 1/2 and b = 1/4 for conventional [1/1] Padé expansion.

To derive a finite-difference representation of the downward continuation operator, we use the following approximation of the exponential operator in equation (9), which is second-order accurate and unitary:

The continuation operator, equation (7), is then approximated as

where, m = 2 b + i a δ' and m^* its complex conjugate. We implemented this algorithm for 2D case, accordingly,

We use a second-order finite-difference scheme to approximate the partial derivative on a regular grid. Using the notation , that is given by:

Substituting this finite-difference approximation into equation (12), we obtained after some algebra the linear system

where N_x is the number of grid points along the x-coordinate. Using the Kronecker symbol δ_i,j, the matrices A_ij and B_ij can be expressed as

with

The matrices A_ij and B_ij in equation (15) are tridiagonal. Solving this linear system we can obtain , the wavefield at z_n + Δz, from the known wavefield at depth z_n. In consequence of equation (9), this scheme can be applied recursively, with appropriate coefficients a_n and b_n, to include higher-order terms in the extrapolation.

Poststack migration

The inputs to the proposed migration method using the complex Padé approximation are the zero-offset data and a reasonable estimate of the velocity distribution. The wavefield is extrapolated in depth, and the imaging condition is accomplished by the exploding-reflector model. The migrated image is obtained at each depth level by adding all frequencies. The output is a migrated depth section in depth domain.

Migration example - SEG/EAGE salt model

We used an exploding-reflector dataset for a 2D slice of the SEG/EAGE salt model to test the migration imaging capability of the complex Padé FD method. The velocity model is shown in Figure 1. The dataset consists of 1290 traces; each trace has 626 time samples with a sampling interval of 8 ms. The velocity model contains 1290 by 300 nodes with a node spacing of 12.192 m along both the horizontal and vertical directions. For comparison, we also apply the split-step Fourier (SSF), Fourier finite-difference (FFD) and phase shift plus interpolation (PSPI) methods to the same dataset. Figures 2 and 3 show images with θ = 0 (real Padé approximation) and 5 degrees, respectively. Some migration artifact is apparent in the result presented in Figure 2. Taking θ ≠ 0 (Fig. 3) we obtain a image with fewer evanescent-energy artifacts than those shown in Figure 2. Figure 4 shows the image obtained with SSF method using the average velocity as a reference velocity. Figure 5 was obtained with PSPI using 5 reference velocities and the Figure 6 shows the result with FFD, now using the minimum velocity as a reference velocity. Comparing the SSF (Fig. 4), PSPI (Fig. 5) and FFD (Fig. 6) and the proposed method (Fig. 3), we see that our method also gives good images of the top and base interfaces of the salt body. Moreover, the image in Figure 3 has fewer numerical artifacts than those obtained with SSF and FFD migration methods.

Prestack migration

Claerbout (1970, 1971) suggested that prestack migration could be done by summing migrated common shot gathers (profiles). Stolt & Benson (1986) applied the method to marine profiles. In order to migrate shot profiles we need the temporal Fourier transform of the recorded common shot profiles at the surface P_U(x, z = 0, ω), the interval velocity field c(x, z) and the shot waveform P_D(x, z = 0, ω).

The method can be summarized by the following steps (Claerbout, 1970, 1971):

1) Downward continuation of the source waveform, P_D:

The downward extrapolation operator,

_s, comes from equation (12) and the downward continuation of the source constitutes a forward modeling.

2) Downward continuation of the waveform recorded by receivers at the surface P_U:

The downward continuation operator for the receivers,

_r, comes from equation (12).

3) Evaluation of the reflectivity function R, in frequency domain, can be formulated by the following equations:

Multiplying both sides of equation (18) by the complex conjugate of P_D, we have:

If equation (19) is inverse Fourier transform over frequency and imaging condition of t = 0 at each z step is invoked, it becomes (Claerbout, 1971):

This imaging principle was proposed by Claerbout (1971) and it can find in space those points where upgoing and downgoing waves are time-coincident, which is where reflectors exist. In order to cure a zero-divide problem, a small positive number Є(z) can be added to the denominator in equation (20), but it is typically difficult to determine since it is data dependent.

The depth image is often obtained by cross-correlating these two wavefields. For several shots it is given by:

2D Marmousi model

In this section we use the 2D Marmousi dataset to test our 2D prestack shot-record depth migration scheme developed based on the complex Padé expansion. The imaging condition is the zero lag of cross-correlation of the downward continued upgoing and downgoing wavefields. The velocity model is illustrated in the Figure 7. Figure 8 shows the migrated sections using the first-order FD method with rotation angles θ = 0º and θ = 10º using only the first term of the equation (8). Figures 9a and 9b further show the migration sections with the second- and third-order FD method for θ = 10º, respectively. The imaging result with the complex Padé propagations shows a significant improvement when compared with the FD imaging results using the real Padé propagations. For comparison we also migrated this dataset with the split-step Fourier (SSF) method using the average velocity as a reference velocity (Fig. 10a) and the phase shift plus interpolation (PSPI) algorithm (Fig. 10b) implemented with 8 reference velocities for each depth extrapolation.

The results of FD migration using complex Padé approximation in Figure 9 are similar to the image SSF and PSPI. Some differences between the images are concentrated upper left part of the image, where the FD image lack focusing for the steep reflectors associated to the faults. We believe these differences can be reduced with a better choice of the rotation angle θ. Further study is still needed to determine an optimum choice of this parameter for each velocity model.

CONCLUSIONS

We have developed a complex Padé FD technique for seismic migration. The method uses the acoustic wave equation directly. Using the Padé approximation of an exponential operator, our method allows the use of a larger grid spacing in depth. A complex treatment of the propagation operator is used to treat the evanescent waves and improve the stability of the numerical approximation. The synthetic results show that high quality images can be obtained with this proposed approach. This technique could be easily extended to the Fourier finite-difference method and also can be improved by using second- or higher-order approximation terms. Further study is still needed to determine the best rotation angle for each velocity model.

ACKNOWLEDGMENTS

This research has been partially supported by the Research Foundation of State of Bahia (FAPESB) and the National Research Council (CNPq) of Brazil. We also thank CT-PETRO/CNPq-FINEP for provided funding.

Recebido em 26 setembro, 2008 / Aceito em 16 janeiro, 2009

Received on September 26, 2008 / Accepted on January 16, 2009

NOTES ABOUT THE AUTHORS

Reynam da Cruz Pestana received his B.Sc. in Physics in 1984 and his PhD in Geophysics in 1988, both degrees from Federal University of Bahia (UFBA). He was a post-doc at Karlsruhe University from 1988 to 1991 and at the University of Texas at Austin from 1998 to 1999 with grants from CNPq. Presently, he is an Associate Professor at the Department of Physics of the Earth and Environment of the Physics Institute of UFBA. He has experience in Geosciences, with emphasis in Applied Geophysics, working mainly in processing and imaging of seismic and GPR data.

Jacira Cristina Batista de Freitas received her B.Sc. in Physics in 1977 and her PhD in Geophysics in 2004, both degrees from Federal University of Bahia (UFBA). Presently, she is an Associate Professor at the Department of Physics of the Earth and Environment of the Physics Institute of UFBA. She has experience in Geosciences, with emphasis in Applied Geophysics, working mainly in theory of GPR data. Her fields of interest include Ray Theory and Gaussian Beam Method, Potential Methods and Mathematical Methods for Geophysicists.

Jessé Carvalho Costa received his diploma in Physics in 1983 from the Physics Department, Federal University of Pará (UFPA) and a Doctor degree in Geophysics in 1993 from the Geophysics Department at the same University. He was a Summer Student at Schlumberger Cambridge Research in 1991 and 1992. He spent 1994 and 1995 as a post-doc in the Stanford Tomography Project at Stanford University. He held a faculty position the Physics Department at UFPA from 1989 to 2003. Currently he is Associate Professor at the Geophysics Department, UFPA. His fields of interest include seismic anisotropy, traveltime tomography and seismic modeling.

• AMAZONAS D, COSTA JC, PESTANA R & SCHLEICHER J. 2007. Wide angle FD and FFD migration using complex Padé approximations, 69^th Meeting. EAGE Conference & Exhibition. CD-ROM.
• CLAERBOUT JF. 1970. Coarse grid calculations of waves in inhomogeneous media with application to delineation of complicated seismic structures. Geophysics, 35: 406-418.
• CLAERBOUT JF. 1971. Toward a unified theory of reflector mapping. Geophysics, 36: 467-481
• CLAERBOUT JF. 1985. Imaging the Earth's interior. Blackwell Scientific Publications.
• GAZDAG J & SGUAZZERO P. 1984. Migration of seismic data by phase shift plus interpolation. Geophysics, 49: 124-131.
• HALE D. 1991. Stable explicit depth extrapolation of seismic wavefields. Geophysics, 56: 1770-1777.
• MILLINAZZO FA, ZALA CA & BROOKE GH. 1997. Square-root approximations for parabolic equation algorithms. J. Acoust. Soc. Am., 101(2): 760-766.
• KOSLOFF D, TAL-EZER H, BARTANA A, RAGOZA E & SHABELANSKY A. 2006. Three dimensional wave equation depth migration by a direction solution method. Expanded Abstracts 76^th Ann. Inter. Mtg., SEG, 2490-2493.
• RISTOW D & RÜHL T. 1994. Fourier finite-difference migration. Geophysics, 59(12): 1882-1893.
• STOFFA PL, FOKKEMA JT, LUNA FREIRE RM & KESSIGER WP. 1990. Split-step Fourier migration. Geophysics, 55: 410-421.
• STOLT RH & BENSON AK. 1986. Seismic Migration: Theory and practice. Geophysical Press.
• ZHANG L, RECTOR JW & HOVERSTEN GM. 2003. Split-step complex Padé migration. Journal of Seismic Exploration, 12: 229-236.
• ZHANG L, RECTOR JW, HOVERSTEN GM & FOMEL S. 2004. Split-step complex Padé-Fourier depth migration. Expanded Abstracts 74^th Ann. Inter. Mtg., SEG. 989-992.

Publication Dates

• Publication in this collection
  14 Aug 2009
• Date of issue
  Mar 2009

History

• Received
  26 Sept 2008
• Accepted
  16 Jan 2009