# Preliminary Investigation of Wavefield Depth Extrapolation by Two-Way Wave Equations

Authors: Bangyu Wu; Ru-Shan Wu; Jinghuai Gao
Journal: International Journal of Geophysics
 (2012)

---

## Abstract

Most of the wavefield downward continuation migration approaches are relying on one-way wave
equations, which move the seismic energy always in one direction along depth. The one-way
downward continuation migrations only use the primaries for imaging and do not treat secondary
reflections recorded on the surface correctly. In this paper, we investigate wavefield depth
extrapolators based on the full acoustic wave equations, which can propagate wave components to
opposite directions. Several two-way wavefield downward continuation propagators are
numerically tested in this study. Recursively implementing of the depth extrapolator makes it
necessary and important to eliminate the unstable wave modes, that is, evanescent waves. For the
laterally varying velocity media, distinction between the propagating and evanescent wave mode
is less clear. We demonstrate that the spatially localized two-way beamlet propagator is an
effective way to remove the evanescent waves while maintain the propagating mode in laterally
inhomogeneous media.


---

## Body

## 1. Introduction

Downward continuation migration calculates the wavefield at greater depth based on the existing wavefield at the shallower depth. For each frequency component, the wavefield can be downward continued recursively from surface to target depth. These algorithms have the flexibility of migrating the seismic data sequentially in depth and frequency, which leads to substantial reduction of both computational and memory requirements. It is particularly advantageous in a migration-velocity-analysis procedure, since one can analyze one layer for each iteration (layer stripping). It also leads to the definition of another useful family of seismic imaging methods—survey-sinking migration [1].Although many different downward continuation algorithms have been proposed, most of them are based on solving the one-way wave equation. The common ground for those methods is splitting the full wave equation into two one-way wave equations that allow for downward or upward wave propagation separately. The directional splitting of the operator suppresses the up-going propagating waves, thus making it difficult, if not impossible, to use the secondary reflections correctly for imaging [2]. The wavefield downward continuation scheme based on the full acoustic wave equation was first introduced by Kosloff and Baysal [3]. For a background with depth-dependent velocity and zero offset source-receiver configuration, the full wave partial differential equation is changed to a second-order ordinary differential equation, and they solved the equation by the fourth-order Runge-Kutta method. Sandberg and Beylkin [4] extended this method to laterally varying media and migrated the prestack seismic data in the source-receiver survey-sinking style. Sandberg et al. [2] modified the algorithm in Sandberg and Beylkin [4] to take a proper account of secondary reflections recorded at the surface for imaging.The recursively implementation of wavefield depth extrapolators makes operator stability a critical issue. In physical reality, only solutions with harmonical oscillations (propagating wave) or exponentially decay (evanescent wave) with depth are present. However, in numerical algorithms, the artificial evanescent waves, which exponentially grow can also be generated and make the extrapolated wavefield grow out of bounds quickly. Usually, the evanescent waves are eliminated to secure the operator stability. For a laterally uniform medium with depth-dependent velocity, the evanescent waves can be suppressed in the frequency wavenumber domain by using the Fourier transform along space and a simple cutoff filter [5]. For the laterally varying velocity field, the identification of the evanescent wave is less clear. To assure the numerical stability, Kosloff and Baysal suggested using the cutoff filter adjusted to the maximum velocity at each depth [3]. However, such a strategy also discarded certain propagating waves generated by the steeply dipping events, resulting in poor imaging of these structures. Sandberg and Beylkin [4] introduced spectral projectors to remove the evanescent wave modes and leave all propagating modes intact in a variable background.We observe that, in the laterally varying velocity field, the propagating and evanescent waves are in fact precisely defined by the velocity at each spatial point. It is difficult to identify different wave modes in the wavenumber domain by using the global Fourier transform assuming one reference velocity at each depth. It leads us to consider the spatially localized wave propagator with reference velocities for different locations. Beamlet propagator [6] provides an alternative approach for seismic wave propagation and imaging based on local velocity perturbation. Wavelet transform is applied along space to shuffle the wavefield between space and local space-wavenumber domain (beamlet domain). Each beamlet is propagated to next depth with local reference velocity. Beamlet propagators have been developed successfully with the Local Cosine Basis (LCB) [7] for the orthogonal and Gabor-Daubechies Frame (GDF) [8] for the nonorthogonal decomposition. The unwanted waves can be removed by an ideal low-pass filter defined by the local reference velocities, leaving the portion of propagating waves increased.In this study, we first test solving the two-way propagator by the 4th-order Rung-Kutta method, the two-way phase-shift method and the implicit two-way phase-shift method, which avoids splitting of the full wave propagator. In the second tests, we see that the two hyperbolas in the surface wavefield are downward continued by the two-way propagators, with each hyperbola moving in opposite direction in the time space seismogram panel. We derive the two-way beamlet propagator for the local background velocity by using LCB transform and test wave propagation and post stack imaging in the lateral high contrast velocity surroundings.


## 2. Full-Wave Downward Continuation Operator in Frequency Wavenumber Domain

### 2.1. Ordinary Differential Two-Way Wave Equation

In a two-dimensional acoustic medium with constant density, the full way wave equation reads(1)∂2u∂x2+∂2u∂z2=1V2∂2u∂t2,
where x and z are the horizontal and vertical coordinates, u(x,z,t) is the pressure field at time t, and V(x,z) is the acoustic velocity. In a lateral homogenous velocity model, equation (1) can be double transformed in x and t to give
(2)∂2û∂z2=-(ω2V2-Kx2)û(Kx,z,ω),
where Kx is the horizontal wavenumber and ω is the temporal frequency.For the migration along depth axis, the second order equation toz can be written as two coupled one order differential wave equations in the depth variable z:
(3)∂û∂z=q,∂q∂z=-(ω2V2-Kx2)û,
where q represents the vertical pressure gradient. Equation (3) can be solved with standard numerical techniques for ordinary differential equations.The full wave equation is of second order on depth and therefore requires two boundary conditions on the surface to initiate the depth extrapolation: the pressure wavefield and its normal derivative. Usually, only the pressure field is recorded on the surface, the other field must be generated with certain assumptions. We estimate the normal derivative of the wavefieldûz(Kx,z=0,ω) using the same approach as in Kosloff and Baysal [3]. Assuming a constant velocity near the surface, the gradient of the wavefield can be obtained from û(Kx,z=0,ω) as
(4)ûz(Kx,z=0,ω)=ikzû(Kx,z=0,ω),
where kz=ω2/V2-Kx2is the vertical wavenumber. Given ûand ûz at surface, equation (3) can be downward continued recursively to any depth.


### 2.2. Phase-Shift Two-Way Depth Extrapolator

The solution to the second-order wave equation (2) can also be written as
(5)ûz(kx,Δz)=12[(ikzû(kx,0)+ûz(kx,0))eikzΔz+(-ikzû(kx,0)+ûz(kx,0))e-ikzΔz],û(kx,Δz)=12[(û(kx,0)+1ikzûz(kx,0))eikzΔz+(û(kx,0)-1ikzûz(kx,0))e-ikzΔz].
Equation (5) is more easily to be understood as two-way wave equation, and we can pick up the up- and down-going terms without difficulty. When downward continue the up- or down-going wave separately, one of the exponential terms can be set to zero. The wavefield gradient is also not needed to decide the wave propagation direction. Equation (5) then reduces to the phase-shift method of Stolt [9].


### 2.3. Implicit Two-Way Phase-Shift Depth Extrapolation

Equation (5) still explicitly splits the two-way wave equation to up- and down-going waves. The phase-shift terms exp(ikzΔz) and exp(-ikzΔz) propagate the up- and down-going components separately.The exponential and the trigonometric functions have the following relations:(6)cos⁡(kzΔz)=eikzΔz+e-ikzΔz2,sin⁡(kzΔz)=eikzΔz-e-ikzΔz2i.
Substitute the exponential terms by the trigonometric functions in (5), we obtain
(7)ûz(kx,Δz)=û(kx,0)(-kzsin⁡(kzΔz))+ûz(kx,0)cos⁡(kzΔz),û(kx,Δz)=û(kx,0)cos⁡(kzΔz)+ûz(kx,0)sin⁡(kzΔz)kz.
Equation (7) avoids splitting the wave equation into the up- and down-going propagators. In practical implementation, (7) is more stable than (5). If kz approaches to zero, kz in denominator of (5) may cause numerical instability. However, in (7), kz and sin(kzΔz) form a sinc function, which eliminates the singularities caused by kz=0.


## 2.1. Ordinary Differential Two-Way Wave Equation

In a two-dimensional acoustic medium with constant density, the full way wave equation reads(1)∂2u∂x2+∂2u∂z2=1V2∂2u∂t2,
where x and z are the horizontal and vertical coordinates, u(x,z,t) is the pressure field at time t, and V(x,z) is the acoustic velocity. In a lateral homogenous velocity model, equation (1) can be double transformed in x and t to give
(2)∂2û∂z2=-(ω2V2-Kx2)û(Kx,z,ω),
where Kx is the horizontal wavenumber and ω is the temporal frequency.For the migration along depth axis, the second order equation toz can be written as two coupled one order differential wave equations in the depth variable z:
(3)∂û∂z=q,∂q∂z=-(ω2V2-Kx2)û,
where q represents the vertical pressure gradient. Equation (3) can be solved with standard numerical techniques for ordinary differential equations.The full wave equation is of second order on depth and therefore requires two boundary conditions on the surface to initiate the depth extrapolation: the pressure wavefield and its normal derivative. Usually, only the pressure field is recorded on the surface, the other field must be generated with certain assumptions. We estimate the normal derivative of the wavefieldûz(Kx,z=0,ω) using the same approach as in Kosloff and Baysal [3]. Assuming a constant velocity near the surface, the gradient of the wavefield can be obtained from û(Kx,z=0,ω) as
(4)ûz(Kx,z=0,ω)=ikzû(Kx,z=0,ω),
where kz=ω2/V2-Kx2is the vertical wavenumber. Given ûand ûz at surface, equation (3) can be downward continued recursively to any depth.


## 2.2. Phase-Shift Two-Way Depth Extrapolator

The solution to the second-order wave equation (2) can also be written as
(5)ûz(kx,Δz)=12[(ikzû(kx,0)+ûz(kx,0))eikzΔz+(-ikzû(kx,0)+ûz(kx,0))e-ikzΔz],û(kx,Δz)=12[(û(kx,0)+1ikzûz(kx,0))eikzΔz+(û(kx,0)-1ikzûz(kx,0))e-ikzΔz].
Equation (5) is more easily to be understood as two-way wave equation, and we can pick up the up- and down-going terms without difficulty. When downward continue the up- or down-going wave separately, one of the exponential terms can be set to zero. The wavefield gradient is also not needed to decide the wave propagation direction. Equation (5) then reduces to the phase-shift method of Stolt [9].


## 2.3. Implicit Two-Way Phase-Shift Depth Extrapolation

Equation (5) still explicitly splits the two-way wave equation to up- and down-going waves. The phase-shift terms exp(ikzΔz) and exp(-ikzΔz) propagate the up- and down-going components separately.The exponential and the trigonometric functions have the following relations:(6)cos⁡(kzΔz)=eikzΔz+e-ikzΔz2,sin⁡(kzΔz)=eikzΔz-e-ikzΔz2i.
Substitute the exponential terms by the trigonometric functions in (5), we obtain
(7)ûz(kx,Δz)=û(kx,0)(-kzsin⁡(kzΔz))+ûz(kx,0)cos⁡(kzΔz),û(kx,Δz)=û(kx,0)cos⁡(kzΔz)+ûz(kx,0)sin⁡(kzΔz)kz.
Equation (7) avoids splitting the wave equation into the up- and down-going propagators. In practical implementation, (7) is more stable than (5). If kz approaches to zero, kz in denominator of (5) may cause numerical instability. However, in (7), kz and sin(kzΔz) form a sinc function, which eliminates the singularities caused by kz=0.


## 3. Two-Way Beamlet Propagators for Local Background Reference Velocities

If the velocity is a function only of depth, the migration in frequency wavenumber domain is accurate and efficient since the Fast Fourier Transform can be used. For a laterally varying medium, the selection of reference velocities is required for the phase-shift operator. To improve the accuracy, part of the computation needs to be performed in the space domain. These approaches can be categorized as mixed domain or dual domain methods. Beamlet migration is a one-way wave equation-based dual domain imaging method. The wave fields are spatially partitioned with local windows and propagated with beamlet propagators (in beamlet domain) followed by local perturbation corrections (in space domain). In local phase space, the evanescent and propagating modes are defined by the local reference velocities and can be processed separately.


### 3.1. Local Cosine Basis and Wavefield Decomposition

The seismic wavefield and its normal derivative can be represented by local Fourier frame, such as the Gabor-Daubechies frame, or the local trigonometric basis (local cosine/sine basis). The local cosine bases are localized version of cosine bases, using overlapped bell functions to decompose signals in separable blocks by fast algorithm [10]. The basis element is characterized by position x̅n, the nominal length of the window Ln=x̅n+1-x̅n, and local wavenumber indexm, (m=0⋯M-1,M denotes the total sample points of the interval) as follows:
(8)bmn(x)=2LnBn(x)cos⁡(π(m+12)x-x̅nLn),
where Bn(x) is a bell function, which is smooth and supported in the compact interval [x̅n-ε,x̅n+1+ε′] for x̅n-ε≤x̅n+1+ε′, where ε, ε′ are the left and right overlapping radius.Wavefieldû(x) and its derivativeûz(x) at depth z can be decomposed into beamlets with windows along the x-axis
(9)ûz(x)=∑n∑m〈ûz(x),bmn(x)〉bmn(x)=∑n∑mUz(x̅n,ξ̅m)bmn(x),û(x)=∑n∑m〈û(x),bmn(x)〉bmn(x)=∑n∑mU(x̅n,ξ̅m)bmn(x).
Here ξ̅m=π(m+1/2)/Ln, U(x̅n,ξ̅m) and Uz(x̅n,ξ̅m) are beamlet coefficients. 〈,〉 stands for inner product.Local cosine bases are real functions. In case of complex wavefield, LCB decomposition is applied to both the real and imaginary parts of the wavefield [7].


### 3.2. LCB Beamlet Domain Two-Way Downward Propagator

For easy implementation, we derive the two-way beamlet propagator starting from (5); we have
(10)ûz(kx,Δz)=12[ikzeikzΔzû(kx,0)+eikzΔzûz(kx,0)+(-ikze-ikzΔz)û(kx,0)+e-ikzΔzûz(kx,0)]=12[P10+û(kx,0)+P00+ûz(kx,0)+P10-û(kx,0)+P00-ûz(kx,0)P10+û],û(kx,Δz)=12[eikzΔzû(kx,0)+eikzΔzikzûz(kx,0)+e-ikzΔzû(kx,0)+e-ikzΔz-ikzûz(kx,0)]=12[P11+û(kx,0)+P01+ûz(kx,0)+P11-û(kx,0)+P01-ûz(kx,0)].
Here Pwith different sub- and superscripts are the phase-shift operator to extrapolate the wavefields and their derivative from depth z to z+Δz. Next, we derive the wave propagator in the LCB beamlet domain following [7]. We define
(11)amn(x)=12π∫dξeiξxPbmn(ξ),
as the evolved LCB beamlet by background propagation using the reference velocity V0(x̅n,z) in window Bn(x). bmn(ξ) is the wavenumber domain beamlet basis. For each exponential term in (10), the evanescent waves are discarded based on the following judgment:
(12)exp⁡(±i(ωV0(x̅n,z))2-kx2Δz)={exp⁡(±i(ωV0(x̅n,z))2-kx2Δz),kx2<(ωV0(x̅n,z))2,0,kx2≥(ωV0(x̅n,z))2.
The redecomposition of the propagated beamlets amn(x) into new beamlets formulates the beamlet propagator, written as
(13)P(ξ̅j,x̅l;ξ̅m,x̅n)=〈amn(x),bjl(ξ)〉.
Thus we obtain the exact expression of the two-way beamlet background propagator components
(14)P10+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)ikzeikzΔz,P00+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)eikzΔz,P10-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)(-ikze-ikzΔz),P00-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)e-ikzΔz,P11+(ξ̅j,x̅l;ξ̅m,x̅n)=P00+(ξ̅j,x̅l;ξ̅m,x̅n),P01+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)eikzΔzikz,P11-(ξ̅j,x̅l;ξ̅m,x̅n)=P00-(ξ̅j,x̅l;ξ̅m,x̅n),P01-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)e-ikzΔz-ikz.
Here the propagator elements P(ξ̅j,x̅l;ξ̅m,x̅n) represent the coupling factors from the beamlet coefficient (ξ̅m,x̅n) to beamlet (ξ̅j,x̅l) at new depth level. Due to the propagation of up- and down-going waves and their derivatives, the two-way beamlet propagator elements are eight times more than that of the one-way propagator. In this algorithm, the computation of two-way wavefield downward extrapolation is approximately 8 times more intensive than the one-way beamlet methods.


## 3.1. Local Cosine Basis and Wavefield Decomposition

The seismic wavefield and its normal derivative can be represented by local Fourier frame, such as the Gabor-Daubechies frame, or the local trigonometric basis (local cosine/sine basis). The local cosine bases are localized version of cosine bases, using overlapped bell functions to decompose signals in separable blocks by fast algorithm [10]. The basis element is characterized by position x̅n, the nominal length of the window Ln=x̅n+1-x̅n, and local wavenumber indexm, (m=0⋯M-1,M denotes the total sample points of the interval) as follows:
(8)bmn(x)=2LnBn(x)cos⁡(π(m+12)x-x̅nLn),
where Bn(x) is a bell function, which is smooth and supported in the compact interval [x̅n-ε,x̅n+1+ε′] for x̅n-ε≤x̅n+1+ε′, where ε, ε′ are the left and right overlapping radius.Wavefieldû(x) and its derivativeûz(x) at depth z can be decomposed into beamlets with windows along the x-axis
(9)ûz(x)=∑n∑m〈ûz(x),bmn(x)〉bmn(x)=∑n∑mUz(x̅n,ξ̅m)bmn(x),û(x)=∑n∑m〈û(x),bmn(x)〉bmn(x)=∑n∑mU(x̅n,ξ̅m)bmn(x).
Here ξ̅m=π(m+1/2)/Ln, U(x̅n,ξ̅m) and Uz(x̅n,ξ̅m) are beamlet coefficients. 〈,〉 stands for inner product.Local cosine bases are real functions. In case of complex wavefield, LCB decomposition is applied to both the real and imaginary parts of the wavefield [7].


## 3.2. LCB Beamlet Domain Two-Way Downward Propagator

For easy implementation, we derive the two-way beamlet propagator starting from (5); we have
(10)ûz(kx,Δz)=12[ikzeikzΔzû(kx,0)+eikzΔzûz(kx,0)+(-ikze-ikzΔz)û(kx,0)+e-ikzΔzûz(kx,0)]=12[P10+û(kx,0)+P00+ûz(kx,0)+P10-û(kx,0)+P00-ûz(kx,0)P10+û],û(kx,Δz)=12[eikzΔzû(kx,0)+eikzΔzikzûz(kx,0)+e-ikzΔzû(kx,0)+e-ikzΔz-ikzûz(kx,0)]=12[P11+û(kx,0)+P01+ûz(kx,0)+P11-û(kx,0)+P01-ûz(kx,0)].
Here Pwith different sub- and superscripts are the phase-shift operator to extrapolate the wavefields and their derivative from depth z to z+Δz. Next, we derive the wave propagator in the LCB beamlet domain following [7]. We define
(11)amn(x)=12π∫dξeiξxPbmn(ξ),
as the evolved LCB beamlet by background propagation using the reference velocity V0(x̅n,z) in window Bn(x). bmn(ξ) is the wavenumber domain beamlet basis. For each exponential term in (10), the evanescent waves are discarded based on the following judgment:
(12)exp⁡(±i(ωV0(x̅n,z))2-kx2Δz)={exp⁡(±i(ωV0(x̅n,z))2-kx2Δz),kx2<(ωV0(x̅n,z))2,0,kx2≥(ωV0(x̅n,z))2.
The redecomposition of the propagated beamlets amn(x) into new beamlets formulates the beamlet propagator, written as
(13)P(ξ̅j,x̅l;ξ̅m,x̅n)=〈amn(x),bjl(ξ)〉.
Thus we obtain the exact expression of the two-way beamlet background propagator components
(14)P10+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)ikzeikzΔz,P00+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)eikzΔz,P10-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)(-ikze-ikzΔz),P00-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)e-ikzΔz,P11+(ξ̅j,x̅l;ξ̅m,x̅n)=P00+(ξ̅j,x̅l;ξ̅m,x̅n),P01+(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)eikzΔzikz,P11-(ξ̅j,x̅l;ξ̅m,x̅n)=P00-(ξ̅j,x̅l;ξ̅m,x̅n),P01-(ξ̅j,x̅l;ξ̅m,x̅n)=12π∫dξbjl*(ξ)bmn(ξ)e-ikzΔz-ikz.
Here the propagator elements P(ξ̅j,x̅l;ξ̅m,x̅n) represent the coupling factors from the beamlet coefficient (ξ̅m,x̅n) to beamlet (ξ̅j,x̅l) at new depth level. Due to the propagation of up- and down-going waves and their derivatives, the two-way beamlet propagator elements are eight times more than that of the one-way propagator. In this algorithm, the computation of two-way wavefield downward extrapolation is approximately 8 times more intensive than the one-way beamlet methods.


## 4. Numerical Experiments

### 4.1. Impulse Responses in Homogenous Media

The impulse responses are computed by different depth extrapolation propagators mentioned above. Considering a domain of 3.75 km deep and 12 km wide, with a constant velocity of 2 km/s, we position a point Ricker wavelet source with dominat-frequency 15 Hz atx=6.4 km and z=0 km. The space and depth sampling intervals are 25 m.The impulse responses by different depth propagators are shown in Figure1 with Figure 1(a) from (3) by the 4th-order Runge-Kutta method, Figure 1(b) two-way phase-shift method (5) and Figure 1(c) the implicit two-way phase-shift method (7), respectively. Figure 1(d) is from the one-way phase-shift method [9]. All the responses look much the same. From these tests, we see that the two-way depth downward continuation methods have the same accuracy as the phase-shift method in the homogenous media as expected.Impulse responses for a homogeneous medium with velocity of 2 km/s from source Ricker wavelet with dominant frequency 15 Hz. (a) Rung-Kutta method; (b) two-way phase-shift method (5); (c) implicit two-way phase-shift depth extrapolation by (7); (d) one-way phase-shift propagator eikzΔz.
(a)(b)(c)(d)During migration, whenω2/V2-Kx2<0, all the wave fields extrapolating by the two-way wave equation blow up rapidly, making it necessary to suppress these unstable modes (evanescent waves) in downward extrapolation algorithms. For a laterally uniform medium, Kosloff and Baysal [3] proposed suppressing the evanescent waves in the wavenumber domain by using a simple cutoff filter (same as used here).


### 4.2. Impulse Responses in Two-Layer Velocity Models

Shown in Figure2 are impulse responses at three time moments for the two-layer models using the implicit two-way phase-shift extrapolator. The velocity change interfaces are in the middle of the depth axis. Left column is for the −50% contrast medium (from 2 km/s to 1 km/s), and the right is for the +50% contrast medium (from 2 km/s to 3 km/s).Impulse responses of the two-way phase shift extrapolator in (7) in two-layer models. Left column is for the −50% contrast velocity model (from 2 km/s to 1 km/s); right column is for the +50% contrast velocity model (from 2 km/s to 3 km/s). Rows from top to bottom are for successive of time sequences.
(a)(b)(c)(d)(e)(f)In Figure2, panels from top to bottom are in order of time sequence. We see that, with the presence of the vertical velocity contrast, there are some anticausal waves enter from the bottoms of the model. In the +50% velocity contrast case, there are vertical artifacts propagating horizontally in the high-velocity region.Figure3 shows the seismograms at depths above and below the interface. The left column is for the −50% velocity contrast and the right column for the +50%. The top panels are seismograms before passing through the interface and the bottom panels, below the interface. We see the hyperbola splits into two after passing the interface and one of which propagates in the opposite direction in the time space panel.Seismograms at depths before (a, b) and after (c, d) pass through the −50% (a, c) and +50% (b, d) contrast velocity interface.
(a)(b)(c)(d)Kosloff and Baysal described the appearance of the up-going energy as the inherent nonuniqueness in the conceptual model on which the migration was based. With the exploding reflector model, a surface recording alone cannot determine the amount of down-going energy at zero time. To do this, a set of geophones need to be placed beneath the structure [3]. To our understanding, the anticasual event appeared after passing through the interfaces because the mathematical downward continuation algorithms do not obey the nature of the wave propagation. The reflected waves are generated when the downward continuation procedure encounters the velocity interface. Physically, the reflected wave should go upward; however, the downward continuation procedure carries it to the future depth by the two-way extrapolator. On the snapshots in Figure 2, the two-way propagators can propagate the up-going wave reflecting by the interface but put it in the wrong spatial location.Seismic imaging based on the two-way depth migration, when high velocity contrast exists, the migrated image will contain low-frequency artifacts similar to the RTM-migrated data. Techniques like smoothing the migration velocity model [2, 4] or applying a Laplacian filter to the final stacked image can effectively dampen the artifacts.


### 4.3. Propagating Wavefield Components in Two Directions in Seismogram

In this section, we test the downward continuation algorithms by the two-way extrapolator with data containing both down-going and up-going waves. Figure4(a) shows a wavefield on the surface containing two hyperbolas, tagged with A and B. The gradient for each hyperbola is calculated separately and added together with opposite signs at the surface, shown in Figure 4(b) (Note the polarity difference between the two gradient fields). Figures 4(c) and 4(d) are seismograms for two depths during downward extrapolation. We see hyperbola A and B moving in opposite directions during depth extrapolation by the two-way wave equations. After some depth, hyperbola A focused into a point, while hyperbola B is continuing to spread. There is no artifacts or cross coupling between the two events. This shows that the two-way extrapolator works in the case of homogeneous or smoothly varying media. For media with sharp boundaries, we can perform some smoothing operation to the media before depth extrapolation or do more study to reduce the artifacts when passing the boundary.Seismograms with two hyperbolas on the surface (a) and its derivative (b); (c) seismogram downward continued to depth 0.55 km; (d) seismogram downward continued to depth 1 km.
(a)(b)(c)(d)


### 4.4. Two-Way Beamlet Propagation and Imaging in Inhomogeneous Media

The beamlet space localization gives the possibility and flexibility in adapting to the fast varying lateral heterogeneities for wave propagation. Theoretically, when there are no lateral velocity variations, the space localized propagators will produce the same result as the methods with only one reference velocity at each depth. However, because the beamlet propagator is the integral of the fine sampling horizontal wavenumbers (14), it is more stable at the vicinity of transiting from propagating modes to evanescent modes. Demonstrating on the snapshots of the impulse responses for the two-layer model, there are less horizontally propagating vertical events in the high-velocity region (shown in the right column of Figure 5).Impulse responses of the two-way beamlet propagator in two-layer models. Left column is for the −50% contrast velocity model (from 2 km/s to 1 km/s); right column is for the +50% contrast velocity model (from 2 km/s to 3 km/s). Rows from top to bottom are for successive of time sequences.
(a)(b)(c)(d)(e)(f)Figure6 shows the background reference velocity for the localized beamlet propagators for 2D SEG/EAGE salt model. The local background velocity reflects the velocity homogeneity in a great degree, such as the salt body and some reflectors. Wave propagation in the local background velocity keeps more propagating waves with high accuracy comparing with the methods only selecting one reference velocity for each depth.Figure 6
Local reference velocities for beamlet depth extrapolation.Figure7 shows impulse responses based on one-way and two-way beamlet propagators in the piecewise local background velocity shown in Figure 6. Panels in left column are for the one-way beamlet propagator and the right column, the two way beamlet propagator. Panels in the same row are at the same time moments for comparison.Impulse responses in a piecewise homogenous media with left and right columns for the one-way and two-way beamlet propagators. From top to bottom are snapshots at times 0.456 s, 0.912 s, and 1.368 s.
(a)(b)(c)(d)(e)(f)The snapshots demonstrate that the localized propagators can handle multipath arrivals with smooth wavefronts amplitude induced by complex velocity functions. We also see that the sharp velocity discontinuity introduce extra signal in the two-way beamlet propagation.In the left column of Figure8, we present snapshots of the two-way beamlet propagator in a smoothed version of the SEG velocity model. Examination of the results in Figures 7(b),  7(d), and 7(f) and Figures 8(a),  8(b), and 8(c) shows that the extra signals produced by the sharp velocity boundaries in the model are greatly suppressed. In Figures 8(b), 8(d), and 8(f), snapshots generated by finite difference method on the velocity in Figure 6 are displayed for comparison.(a, c, e) are snapshots of the two-way beamlet propagator on the smoothed version of the SEG velocity model; (b, d, f) are snapshots on velocity model shown in Figure6 by finite difference method.
(a)(b)(c)(d)(e)(f)Figure9 illustrates the post stack imaging performance of the localized propagator on 2D SEG/EAGE synthetic data sets. Figure 9(a) shows the original velocity model, with minimum and maximum velocities are 1500 m/s and 4410 m/s. Figure 9(b) shows the post stack image of the one-way beamlet migration with local perturbation correction (for details of local perturbation correction see reference [7]).2D SEG/EAGE post stack imaging by using different propagators. (a) Original velocity model; (b) one-way beamlet migration with local perturbation correction; (c) one-way beamlet migration with local reference velocity shown in Figure6; (d) two-way beamlet migration with local reference velocity shown in Figure 6.
(a)(b)(c)(d)Figures9(c) and 9(d) display the image based on one-way and two-way beamlet propagators using the local reference velocities, without perturbation correction. In Figure 9(c), the boundary of the salt body and most of the reflectors are all imaged. Due to the accumulated errors from velocity perturbations, the image is unsatisfactory in the sub-salt area. In presence of strong velocity contrasts, the anticausal artifacts generated by the two-way extrapolator will significantly contaminate the final image. Kosloff and Baysal eliminated the anti-causal energy from the depth section by filtering out components with negative vertical wavenumbers [3].A blurred version of the original velocity model can also be used as the background velocity to avoid the migration artifacts. Figure10(a) is a smoothed SEG velocity and Figure 10(b) is corresponding local background velocity. Image in Figure 10(c) shows that some of the artifacts in Figure 9(d) are removed. For the smoothed velocity, because the selections constant reference velocities in each window, the background velocity for the two-way beamlet propagators is still piecewise and there are still some remaining artifacts in the migrated image.2D SEG/EAGE post stack imaging on the smoothed version of the velocity model. (a) Velocity model obtained by blurring the original model; (b) local background velocity for migration by the two-way beamlet propagator; (c) migrated image using the background velocity in (b).
(a)(b)(c)


## 4.1. Impulse Responses in Homogenous Media

The impulse responses are computed by different depth extrapolation propagators mentioned above. Considering a domain of 3.75 km deep and 12 km wide, with a constant velocity of 2 km/s, we position a point Ricker wavelet source with dominat-frequency 15 Hz atx=6.4 km and z=0 km. The space and depth sampling intervals are 25 m.The impulse responses by different depth propagators are shown in Figure1 with Figure 1(a) from (3) by the 4th-order Runge-Kutta method, Figure 1(b) two-way phase-shift method (5) and Figure 1(c) the implicit two-way phase-shift method (7), respectively. Figure 1(d) is from the one-way phase-shift method [9]. All the responses look much the same. From these tests, we see that the two-way depth downward continuation methods have the same accuracy as the phase-shift method in the homogenous media as expected.Impulse responses for a homogeneous medium with velocity of 2 km/s from source Ricker wavelet with dominant frequency 15 Hz. (a) Rung-Kutta method; (b) two-way phase-shift method (5); (c) implicit two-way phase-shift depth extrapolation by (7); (d) one-way phase-shift propagator eikzΔz.
(a)(b)(c)(d)During migration, whenω2/V2-Kx2<0, all the wave fields extrapolating by the two-way wave equation blow up rapidly, making it necessary to suppress these unstable modes (evanescent waves) in downward extrapolation algorithms. For a laterally uniform medium, Kosloff and Baysal [3] proposed suppressing the evanescent waves in the wavenumber domain by using a simple cutoff filter (same as used here).


## 4.2. Impulse Responses in Two-Layer Velocity Models

Shown in Figure2 are impulse responses at three time moments for the two-layer models using the implicit two-way phase-shift extrapolator. The velocity change interfaces are in the middle of the depth axis. Left column is for the −50% contrast medium (from 2 km/s to 1 km/s), and the right is for the +50% contrast medium (from 2 km/s to 3 km/s).Impulse responses of the two-way phase shift extrapolator in (7) in two-layer models. Left column is for the −50% contrast velocity model (from 2 km/s to 1 km/s); right column is for the +50% contrast velocity model (from 2 km/s to 3 km/s). Rows from top to bottom are for successive of time sequences.
(a)(b)(c)(d)(e)(f)In Figure2, panels from top to bottom are in order of time sequence. We see that, with the presence of the vertical velocity contrast, there are some anticausal waves enter from the bottoms of the model. In the +50% velocity contrast case, there are vertical artifacts propagating horizontally in the high-velocity region.Figure3 shows the seismograms at depths above and below the interface. The left column is for the −50% velocity contrast and the right column for the +50%. The top panels are seismograms before passing through the interface and the bottom panels, below the interface. We see the hyperbola splits into two after passing the interface and one of which propagates in the opposite direction in the time space panel.Seismograms at depths before (a, b) and after (c, d) pass through the −50% (a, c) and +50% (b, d) contrast velocity interface.
(a)(b)(c)(d)Kosloff and Baysal described the appearance of the up-going energy as the inherent nonuniqueness in the conceptual model on which the migration was based. With the exploding reflector model, a surface recording alone cannot determine the amount of down-going energy at zero time. To do this, a set of geophones need to be placed beneath the structure [3]. To our understanding, the anticasual event appeared after passing through the interfaces because the mathematical downward continuation algorithms do not obey the nature of the wave propagation. The reflected waves are generated when the downward continuation procedure encounters the velocity interface. Physically, the reflected wave should go upward; however, the downward continuation procedure carries it to the future depth by the two-way extrapolator. On the snapshots in Figure 2, the two-way propagators can propagate the up-going wave reflecting by the interface but put it in the wrong spatial location.Seismic imaging based on the two-way depth migration, when high velocity contrast exists, the migrated image will contain low-frequency artifacts similar to the RTM-migrated data. Techniques like smoothing the migration velocity model [2, 4] or applying a Laplacian filter to the final stacked image can effectively dampen the artifacts.


## 4.3. Propagating Wavefield Components in Two Directions in Seismogram

In this section, we test the downward continuation algorithms by the two-way extrapolator with data containing both down-going and up-going waves. Figure4(a) shows a wavefield on the surface containing two hyperbolas, tagged with A and B. The gradient for each hyperbola is calculated separately and added together with opposite signs at the surface, shown in Figure 4(b) (Note the polarity difference between the two gradient fields). Figures 4(c) and 4(d) are seismograms for two depths during downward extrapolation. We see hyperbola A and B moving in opposite directions during depth extrapolation by the two-way wave equations. After some depth, hyperbola A focused into a point, while hyperbola B is continuing to spread. There is no artifacts or cross coupling between the two events. This shows that the two-way extrapolator works in the case of homogeneous or smoothly varying media. For media with sharp boundaries, we can perform some smoothing operation to the media before depth extrapolation or do more study to reduce the artifacts when passing the boundary.Seismograms with two hyperbolas on the surface (a) and its derivative (b); (c) seismogram downward continued to depth 0.55 km; (d) seismogram downward continued to depth 1 km.
(a)(b)(c)(d)


## 4.4. Two-Way Beamlet Propagation and Imaging in Inhomogeneous Media

The beamlet space localization gives the possibility and flexibility in adapting to the fast varying lateral heterogeneities for wave propagation. Theoretically, when there are no lateral velocity variations, the space localized propagators will produce the same result as the methods with only one reference velocity at each depth. However, because the beamlet propagator is the integral of the fine sampling horizontal wavenumbers (14), it is more stable at the vicinity of transiting from propagating modes to evanescent modes. Demonstrating on the snapshots of the impulse responses for the two-layer model, there are less horizontally propagating vertical events in the high-velocity region (shown in the right column of Figure 5).Impulse responses of the two-way beamlet propagator in two-layer models. Left column is for the −50% contrast velocity model (from 2 km/s to 1 km/s); right column is for the +50% contrast velocity model (from 2 km/s to 3 km/s). Rows from top to bottom are for successive of time sequences.
(a)(b)(c)(d)(e)(f)Figure6 shows the background reference velocity for the localized beamlet propagators for 2D SEG/EAGE salt model. The local background velocity reflects the velocity homogeneity in a great degree, such as the salt body and some reflectors. Wave propagation in the local background velocity keeps more propagating waves with high accuracy comparing with the methods only selecting one reference velocity for each depth.Figure 6
Local reference velocities for beamlet depth extrapolation.Figure7 shows impulse responses based on one-way and two-way beamlet propagators in the piecewise local background velocity shown in Figure 6. Panels in left column are for the one-way beamlet propagator and the right column, the two way beamlet propagator. Panels in the same row are at the same time moments for comparison.Impulse responses in a piecewise homogenous media with left and right columns for the one-way and two-way beamlet propagators. From top to bottom are snapshots at times 0.456 s, 0.912 s, and 1.368 s.
(a)(b)(c)(d)(e)(f)The snapshots demonstrate that the localized propagators can handle multipath arrivals with smooth wavefronts amplitude induced by complex velocity functions. We also see that the sharp velocity discontinuity introduce extra signal in the two-way beamlet propagation.In the left column of Figure8, we present snapshots of the two-way beamlet propagator in a smoothed version of the SEG velocity model. Examination of the results in Figures 7(b),  7(d), and 7(f) and Figures 8(a),  8(b), and 8(c) shows that the extra signals produced by the sharp velocity boundaries in the model are greatly suppressed. In Figures 8(b), 8(d), and 8(f), snapshots generated by finite difference method on the velocity in Figure 6 are displayed for comparison.(a, c, e) are snapshots of the two-way beamlet propagator on the smoothed version of the SEG velocity model; (b, d, f) are snapshots on velocity model shown in Figure6 by finite difference method.
(a)(b)(c)(d)(e)(f)Figure9 illustrates the post stack imaging performance of the localized propagator on 2D SEG/EAGE synthetic data sets. Figure 9(a) shows the original velocity model, with minimum and maximum velocities are 1500 m/s and 4410 m/s. Figure 9(b) shows the post stack image of the one-way beamlet migration with local perturbation correction (for details of local perturbation correction see reference [7]).2D SEG/EAGE post stack imaging by using different propagators. (a) Original velocity model; (b) one-way beamlet migration with local perturbation correction; (c) one-way beamlet migration with local reference velocity shown in Figure6; (d) two-way beamlet migration with local reference velocity shown in Figure 6.
(a)(b)(c)(d)Figures9(c) and 9(d) display the image based on one-way and two-way beamlet propagators using the local reference velocities, without perturbation correction. In Figure 9(c), the boundary of the salt body and most of the reflectors are all imaged. Due to the accumulated errors from velocity perturbations, the image is unsatisfactory in the sub-salt area. In presence of strong velocity contrasts, the anticausal artifacts generated by the two-way extrapolator will significantly contaminate the final image. Kosloff and Baysal eliminated the anti-causal energy from the depth section by filtering out components with negative vertical wavenumbers [3].A blurred version of the original velocity model can also be used as the background velocity to avoid the migration artifacts. Figure10(a) is a smoothed SEG velocity and Figure 10(b) is corresponding local background velocity. Image in Figure 10(c) shows that some of the artifacts in Figure 9(d) are removed. For the smoothed velocity, because the selections constant reference velocities in each window, the background velocity for the two-way beamlet propagators is still piecewise and there are still some remaining artifacts in the migrated image.2D SEG/EAGE post stack imaging on the smoothed version of the velocity model. (a) Velocity model obtained by blurring the original model; (b) local background velocity for migration by the two-way beamlet propagator; (c) migrated image using the background velocity in (b).
(a)(b)(c)


## 5. Discussion and Conclusions

The potential of using multiple reflections for imaging is very attractive for the depth extrapolation based on two-way wave equations. In this paper, we did some preliminary tests for wave propagation and imaging based on the two-way depth extrapolators, which can propagate the wavefield components in opposite directions. The downward continuation scheme is equivalent to an initial value problem and two initial conditions are needed on the surface to start the propagation. The evanescent waves must be eliminated for stable extrapolation. In lateral inhomogeneous media, the identification and elimination of unwanted waves directly relate to the wavefield accuracy and imaging quality. Beamlet propagator extrapolates the wavefield in space-localized phase space. In laterally inhomogeneous media, the evanescent and propagating modes are defined by the local reference velocity in each window, which provides an effective way to improve the accuracy of wave mode identification. The two-way beamlet propagator is computationally more expensive than the one-way propagator. For migrations in the models with sharp boundaries, we must take measures to eliminate anti-casual events generated by the two-way depth propagators passing through the interfaces. More study is needed to develop a depth extrapolation procedure, which can take advantage of the two-way wave propagator.



---

*Source: 968090-2012-06-13.xml*