# Rational approximations for implicit schemes

Series Investigations in Geophysics Öz Yilmaz http://dx.doi.org/10.1190/1.9781560801580 ISBN 978-1-56080-094-1 SEG Online Store

The scalar wave equation is a two-way wave equation in depth that describes propogation of both upcoming and downgoing waves. If we consider the resulting wavefield from the exploding reflectors model as the upcoming waves, then we are really interested in a one-way wave equation to downward continue the upcoming waves. In fact, we normally use some rational approximation to the one-way wave equation in finite-difference implementations.

To get the actual differential equation to be used in downward extrapolation of the upcoming waves, and therefore to perform a finite-difference migration, the general strategy is as follows:

 ${\displaystyle {\frac {\partial ^{2}P}{\partial x^{2}}}+{\frac {\partial ^{2}P}{\partial z^{2}}}-{\frac {1}{v^{2}(x,z)}}{\frac {\partial ^{2}P}{\partial t^{2}}}=0,}$ (12)
where x and z are the space variables, t is the time variable, v is the velocity of wave propagation, and P(x, z, t) is the pressure wavefield.
• Assume constant velocity and perform 3-D Fourier transform of the pressure wavefield. This is equivalent to substituting the plane-wave solution exp(ikxx + ikzziωt) to equation (12). The substitution yields the dispersion relation between the transform variables

 ${\displaystyle k_{z}=\mp {\sqrt {{\frac {\omega ^{2}}{v^{2}}}-k_{x}^{2}}},}$ (13a)
where kx and kz are the wavenumbers in the x and z directions, and ω is the angular temporal frequency.
• We are interested in upcoming waves, hence we only need one of the two solutions. We also want to invoke the exploding reflector model by replacing v by v/2 to obtain the following paraxial dispersion relation

 ${\displaystyle k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}}},}$ (13b)
where the horizontal wavenumber kx has been normalized with respect to 2ω/v.
• Make a rational approximation to the square-root expression in equation (13b) so as to derive a differential equation (Sections D.3 and D.4). This approximation imposes a dip limitation to the differential equation. One approximation to the dispersion relation given by equation (13b) is obtained by Taylor expansion of the square root and retaining the first two terms in the series (Section D.3)

 ${\displaystyle k_{z}={\frac {2\omega }{v}}\left[1-{\frac {1}{2}}\left({\frac {vk_{x}}{2\omega }}\right)^{2}\right].}$ (14a)
This dispersion equation is known as the 15-degree approximation and is the basis for the first finite-difference time migration algorithm developed by Claerbout and Doherty [1]. Albeit no longer in use, we shall review the 15-degree finite-difference algorithm for its historical significance.
• Operate on the pressure wavefield P(kx, kz, ω) with the approximate form of the dispersion relation given by equation (14a), and inverse Fourier transform in the z direction to get the differential form of the approximate one-way wave equation.

 ${\displaystyle {\frac {\partial }{\partial z}}P(k_{x},z,\omega )=-i{\frac {2\omega }{v}}\left[1-{\frac {1}{2}}\left({\frac {vk_{x}}{2\omega }}\right)^{2}\right]P(k_{x},z,\omega ).}$ (14b)
• Recall from Figure 4.1-17 that, after each downward-continuation step, we retard the wave-field by translating it in time so that after migration, events appear in their correct depth locations. The time retardation is done by applying a linear phase shift to the pressure wavefield P

 ${\displaystyle Q=P\exp(-i\omega \tau ),}$ (15a)
where the retarded time is

 ${\displaystyle \tau =\int _{0}^{z}{\frac {dz}{{\bar {v}}(z)}},}$ (15b)
and Q is the retarded wavefield. The velocity ${\displaystyle {\bar {v}}(z)}$ is the horizontally averaged v(x, z). Substitute equation (15a) into (14b) to obtain the differential equation associated with the 15-degree finite-difference algorithm in two parts

 ${\displaystyle {\frac {\partial ^{2}Q}{\partial z\partial t}}={\frac {v}{4}}{\frac {\partial ^{2}Q}{\partial x^{2}}},}$ (16a)

and

 ${\displaystyle {\frac {\partial Q}{\partial z}}=2\left[{\frac {1}{{\bar {v}}(z)}}-{\frac {1}{v(x,z)}}\right]{\frac {\partial Q}{\partial t}},}$ (16b)
where Q is the retarded wavefield. Derivation of equations (16a,16b) is based on the assumption that velocity varies vertically. Nevertheless, in practice, the velocity function in equations (16a,16b) can be varied laterally, provided the variation is smooth. Equation (16a) accounts for collapsing diffraction energy to the apex of the traveltime curve only. Hence, it is referred to as the diffraction term. When lateral velocity variations are significant, the diffraction curve is somewhat like a skewed hyperbola with its apex shifted laterally away from the diffraction source. This lateral shift is accounted for by the thin-lens term given by equation (16b) (Section D.3). If the lateral velocity variations are significant, then the thin-lens term is not negligible. Migration algorithms that implement both the diffraction and thin-lens terms represented by equations (16a,16b) generally are two-step schemes that alternately solve these two terms. To propagate one depth step, first apply the diffraction term on wavefield Q. The thin-lens term then is applied to the output from the diffraction calculation. A migration method that includes the effects of the thin-lens term is called depth migration, since the output section is in depth. Depth migration is warranted if there are strong lateral variations in velocity; in this case, the coefficient of the thin-lens term cannot be negligible. If we assume that velocity varies only in the vertical direction, then ${\displaystyle {\bar {v}}(z)=v(x,\ z).}$ This makes the thin-lens term of equation (16b) vanish, and we are left with the diffraction term of equation (16a). A migration method that implements the diffraction term (equation 16a), only, is known as time migration, the output of which is in time τ of equation (15b). When recast in terms of the τ variable, equation (16a) takes the form

 ${\displaystyle {\frac {\partial ^{2}Q}{\partial \tau \partial t}}={\frac {v^{2}}{8}}{\frac {\partial ^{2}Q}{\partial x^{2}}}.}$ (17)
This is the parabolic equation for time migration.
• Finally, write down the difference forms of the differential operators either in implicit form to be used in finite-difference solution of the parabolic equation (17) for migration.

Boundary and initial conditions are needed to solve the differential equations. The initial condition for migration is the recorded wavefield at the surface z = 0. Also, in migration we assume that the wavefield is zero after a maximum observation time, typically the end time of the recorded trace. Then there are the side boundaries, beyond which assumptions must be made about the form of the wavefield.

In the (x, z, t) coordinates, the seismic section is represented by the x − t plane, while the migrated section (earth) is represented by the x − z plane. Finite-difference migration, as discussed here, extrapolates the x − t plane in finite increments of z and outputs the wavefield at t = 0 at each step (Figure 4.1-19).

There are two ways to downward continue the wavefield recorded at the surface in the computer (Figure 4.1-20). Starting with the wavefield at the surface z = 0 represented by the vectors in x — s1, s2, s3, …, which are perpendicular to the page, we can compute the wavefield at different depth levels using the order of the computation shown in Figure 4.1-20a. Assume zero value for the bottom of the extrapolated wavefield at each depth step. So, for example, using s7, s8, and 0, compute the wavefield at position 1. Then use s6, s7, and the wavefield already computed at position 1, compute that at position 2, and so on. Notice that, in this scheme, we compute the wavefield at all times for one depth step, then compute the wavefield at all times for the next depth step, followed by the next depth step, and so on. Hence, this is called the z-outer computational scheme.

The alternate scheme involves a different order of computation (Figure 4.1-20b). First, compute the wave-field at one time for all depths, then using those already computed values, compute the wavefield at the next shallower time for all depths, and so on. Hence this is called the t-outer computational scheme.

In both schemes, the output migrated section is obtained by collecting the diagonal elements. Depending on the depth step size, which can be conveniently defined as the number of time samples, one collects one or more samples at each depth level. In the example shown in Figure 4.1-20, there are two samples from each depth step collected into the migrated section.

## References

1. Claerbout and Doherty, 1972, Claerbout, J.F. and Doherty, S.M., 1972, Downward continuation of moveout-corrected seismo-grams: Geophysics, 37, 741–768.
2. Claerbout (1976), Claerbout, J.F., 1976, Fundamentals of geophysical data processing: McGraw-Hill Book Co.