Diffraction and ray theory for wave propagation

Seismic Data Analysis

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

H.1 The Kirchhoff integral

Kirchhoff’s integral solution to the scalar wave equation

$\left[{\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial y^{2}}}+{\frac {\partial ^{2}}{\partial z^{2}}}-{\frac {1}{v^{2}(x,y,z)}}{\frac {\partial ^{2}}{\partial t^{2}}}\right]P(x,y,z;t)=0$

(H-1)

is a mathematical statement of Huygen’s principle. In equation (H-1), P(x, y, z; t) is the pressure wavefield propagating in a medium with velocity v(x, y, z). Huygen’s principle states that the pressure disturbance at time t + Δt is the superposition of the spherical waves generated by point sources at time t.

Consider the geometry in Figure H-1 of a point diffractor at a location S(x, y, z) and an observation surface A for the diffraction wavefield generated by the source at the diffractor. Actually, the surface area A is only a portion of a closed surface and is the aperture of the observation made over that closed surface. For convenience, we shall choose the receiver location R(0, 0, 0) on the observation surface A to be at the origin of our coordinate system.

Also for convenience, we shall apply Fourier transform to the wavefield in the time direction

$P(x,y,z,;\omega )=\int P(x,y,z;t)\ \ \exp(-i\omega t)\ dt,$

(H-2a)

where ω is the angular frequency. The inverse transform is given by

$P(x,y,z,;t)=\int P(x,y,z;\omega )\ \ \exp(i\omega t)\ d\omega .$

(H-2b)

Now apply Fourier transform to equation (H-1) in the time direction to obtain

$\left(\nabla ^{2}+{\frac {\omega ^{2}}{v^{2}}}\right)P(x,y,z;\omega )=0,$

(H-3)

where $\nabla ^{2}$ is the Laplacian operator

$\nabla ^{2}P=\left({\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial y^{2}}}+{\frac {\partial ^{2}}{\partial z^{2}}}\right)P.$

We may intuitively state that what we observe at the surface A is what is generated at the source S. This statement is mathematically expressed by Gauss’s divergence theorem as

$\int _{V}\nabla ^{2}PdV=\int _{A}{\frac {\partial P}{\partial n}}\ dA,$

(H-4)

where V is the volume of the region enclosed by the surface A and the derivative ∂P/∂n is taken normal to the surface A in the outward direction. Note that Gauss’s divergence theorem expressed by equation (H-4) transforms a surface integral into a volume integral. For our specific application, this means imaging the earth’s interior from seismic waves observed at the earth’s surface.

We shall solve equation (H-3) for each frequency component ω and sum the resulting solutions over all frequency components to compute the wavefield at the source P(x, y, z; t = 0). The solution obtained by Kirchhoff in 1882 requires Green’s function which describes the propagation outward from a point source with spherical symmetry as

$G(r,\omega )={\frac {1}{r}}\ \ \exp \left(-i{\frac {\omega }{v}}r\right),$

(H-5a)

where

$r={\sqrt {x^{2}+y^{2}+z^{2}}}$

(H-5b)

is the distance between the observation point and the source location. Green’s function given by equation (H-5a) is also a valid solution to equation (H-3):

$\left(\nabla ^{2}+{\frac {\omega ^{2}}{v^{2}}}\right)G(x,y,z;\omega )=0.$

(H-6)

Figure H-1  Geometry of a point diffractor to derive the Kirchhoff integral solution to the scalar wave equation.

We now rewrite equation (H-4) by multiplying both sides with Green’s function G of equation (H-5a) as

$\int _{V}G\nabla ^{2}PdV=\int _{A}G{\frac {\partial P}{\partial n}}\ dA,$

(H-7a)

and rewrite, in turn, equation (H-7a) by interchanging our wave function P with Green’s function G to obtain

$\int _{V}P\nabla ^{2}GdV=\int _{A}P{\frac {\partial G}{\partial n}}dA.$

(H-7b)

Now subtract equation (H-7b) from (H-7a) to get

$\int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=\int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA.$

(H-8)

Equation (H-8) is known as Green’s theorem. Albeit the algebraic steps that we took from Gauss’s theorem given by equation (H-4) to Green’s theorem given by equation (H-8) may appear informal, the intention here is to provide the basic concept underlying the derivation.

Substitute equations (H-3) and (H-6) into the left-hand side of equation (H-8)

$\int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=\int _{V}\left(-G{\frac {\omega ^{2}}{v^{2}}}P+P{\frac {\omega ^{2}}{v^{2}}}G\right)\ dV,$

and hence note that

$\int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=0.$

(H-9)

Now we turn our attention to the right-hand side of equation (H-8). Because Green’s function defined by equation (H-5a) becomes infinite at the source location S, we need to place it inside an infinitesimally small enclosed surface E. This would then require computing the right-hand side of equation (H-8) in two parts — once for the surface E and once for the surface A.

Substitute equation (H-5a) into the right-hand side of equation (H-8) and note, from Figure H-1, that for the surface E, (∂/∂n) = −(∂/∂r):

$\int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=\int _{E}\left[-{\frac {1}{r}}\ \exp {\Bigg (}-i{\frac {\omega }{v}}r{\Bigg )}{\frac {\partial P}{\partial r}}+P{\frac {\partial }{\partial r}}\left({\frac {1}{r}}\ \exp {\Bigg (}-i{\frac {\omega }{v}}r{\Bigg )}\right)\right]\ r^{2}d\Omega ,$

(H-10a)

where dE = r²dΩ and Ω is the solid angle around the point source S in Figure H-1. Carry out the differentiation with respect to r and simplify the right-hand side of equation (H-10a) to obtain

$\int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=-\int _{E}\exp {\bigg (}-i{\frac {\omega }{v}}r{\bigg )}\left(r{\frac {\partial P}{\partial r}}+P+i{\frac {\omega }{v}}rP\right)\ d\Omega .$

(H-10b)

Finally, take the limit r → 0 to obtain the contribution of the surface E:

$\int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=-4\pi P.$

(H-11)

Now substitute, again, equation (H-5a) into the right-hand side of equation (H-8) and note from Figure H-1 that for the surface A, (∂/∂n) = −(∂/∂z):

$\int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA=\int _{A}\left\{-{\frac {1}{r}}\ \exp \left(-i{\frac {\omega }{v}}r\right){\frac {\partial P}{\partial z}}+P{\frac {\partial }{\partial z}}\left[{\frac {1}{r}}\ \exp \left(-i{\frac {\omega }{v}}r\right)\right]\right\}\ dA.$

(H-12)

Carry out the differentiation with respect to z while noting from Figure H-1 that ∂r/∂z = cos θ, and simplify the right-hand side of equation (H-12) to get

$\int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA=-\int _{A}\exp \left(-i{\frac {\omega }{v}}r\right)\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ dA.$

(H-13)

The total contribution to the right-hand side of equation (H-8) actually is the sum of equations (H-11) and (H-13). The left-hand side of equation (H-8) vanishes by way of equation (H-9). Hence, the resulting expression from equation (H-9) is

$4\pi P=\int _{A}\exp \left(-i{\frac {\omega }{v}}r\right)\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ dA.$

(H-14)

Recall that P = P(x, y, z; ω) in equation (H-14) and multiply both sides by exp(iωt) and integrate over frequency ω. The left-hand side of the equation then becomes P(x, y, z; t) by way of the inverse Fourier transform as in equation (H-2b). Thus the resulting expression is

$P(x,y,z;t)={\frac {1}{4\pi }}\int _{\omega }\int _{A}\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ \exp \left[-i\omega \left(t-{\frac {r}{v}}\right)\right]\ dA\ d\omega .$

(H-15)

Define the variable τ = t − r/v as retarded time. Recall from Appendix A that if P(ω) is the Fourier transform of P(t), then exp(−iωr/v)P(ω) is the Fourier transform of P(τ = t − r/v). Also, if P(ω) is the Fourier transform of P(t), then iωP(ω) is the Fourier transform of ∂P/∂t. Incorporate these relations into equation (H-15) and apply inverse Fourier transform to the right-hand side to obtain

$P(x,y,z;\tau )={\frac {1}{4\pi }}\int _{A}\left\{{\frac {1}{r}}\left[{\frac {\partial P}{\partial z}}\right]+{\frac {\cos \theta }{r^{2}}}[P]+{\frac {\cos \theta }{vr}}\left[{\frac {\partial P}{\partial t}}\right]\right\}\ dA,$

(H-16)

where [P] means that the integration over the area A is done using the wavefield P at retarded time τ = t − r/v.

The first term depends on the vertical gradient of the wavefield ∂P/∂z. The second term is called the near-field term since it decays with 1/r². Both terms are neglected in seismic migration. The remaining third term is called the far-field term and it is the foundation of Kirchhoff migration. Writing it in discrete form, we have

$P_{out}={\frac {\Delta x\Delta y}{4\pi }}\sum _{A}{\frac {\cos \theta }{vr}}{\frac {\partial }{\partial t}}P_{in},$

(H-17)

where Δx and Δy are inline and crossline trace spacings, P_out = P(x_out, y_out, z; τ = 2z/v) is the output of migration using the input wavefield P_in = P(x_in, y_in, z = 0; τ = t − r/v) within an areal aperture A.

H.2 The eikonal equation

The eikonal equation is a ray-theoretical approximation to the scalar wave equation. Its solution represents wavefronts of constant phase. The wave is propagated from one wavefront to the next by way of raypaths which are perpendicular to the wavefronts.

Consider a compressional plane wave P(x, y, z; t) in 3-D Cartesian coordinates (x, y, z)

$P(x,y,z;t)=P_{0}\ \exp(-i\omega t+ik_{x}x+ik_{y}y+ik_{z}z),$

(H-18)

where P₀ is the wave amplitude, t is the traveltime, and k_x, k_y, k_z and ω are the Fourier duals of the variables x, y, z and t, respectively. That this equation is a solution to the 3-D scalar wave equation

${\frac {\partial ^{2}P}{\partial x^{2}}}+{\frac {\partial ^{2}P}{\partial y^{2}}}+{\frac {\partial ^{2}P}{\partial z^{2}}}={\frac {1}{v^{2}(x,y,z)}}{\frac {\partial ^{2}P}{\partial t^{2}}}$

(H-19)

can be verified by computing the partial derivatives of the wavefield P(x, y, z; t) given by equation (H-18) and direct substitution into equation (H-19). The result is the dispersion relation

$k_{x}^{2}+k_{y}^{2}+k_{z}^{2}={\frac {\omega ^{2}}{v^{2}}}$

(H-20)

of the scalar wave equation (H-19). Since equation (H-18) satisfies this relation, it is a valid solution to the scalar wave equation. In equation (H-19), v(x, y, z) is the propagation velocity of the compressional plane wave P(x, y, z; t).

We shall rewrite the phase term in the plane-wave solution given by equation (H-18) as

$P(x,y,z;t)=P_{0}\ \exp \left\{-i\omega \left[t-\left({\frac {k_{x}}{\omega }}x+{\frac {k_{y}}{\omega }}y+{\frac {k_{z}}{\omega }}z\right)\right]\right\},$

(H-21)

so that we can define a 3-D traveltime surface T(x, y, z) as

$T(x,y,z)={\frac {k_{x}}{\omega }}x+{\frac {k_{y}}{\omega }}y+{\frac {k_{z}}{\omega }}z.$

(H-22)

Substitute this definition into equation (H-21) to obtain the expression for the plane wave in terms of the traveltime surface T(x, y, z)

$P(x,y,z;t)=P_{0}\ \exp {\Big \{}-i\omega {\Big [}t-T(x,y,z){\Big ]}{\Big \}}.$

(H-23)

Now, we need to verify that this form of the plane wave solution satisfies the scalar wave equation (H-19). Compute the partial derivative of the wave function P(x, y, z; t) of equation (H-23) with respect to the space variable x

${\frac {\partial P}{\partial x}}=P_{0}(i\omega )\left({\frac {\partial T}{\partial x}}\right)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}$

(H-24)

and differentiate the result once more to get

${\frac {\partial ^{2}P}{\partial x^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial x}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial x^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-25a)

Repeat the differentiation for the space variable y

${\frac {\partial ^{2}P}{\partial y^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial y}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial y^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}$

(H-25b)

and for the space variable z

${\frac {\partial ^{2}P}{\partial z^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial z}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial z^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-25c)

Next compute the partial derivative of the wave function P(x, y, z; t) of equation (H-23) with respect to the time variable t

${\frac {\partial ^{2}P}{\partial t^{2}}}=-P_{0}\omega ^{2}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-25d)

Substitute equations (H-25a,H-25b,H-25c,H-25d) into equation (H-19) and combine the terms into real and imaginary parts

$\omega ^{2}\left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-i\omega \left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)={\frac {\omega ^{2}}{v^{2}(x,y,z)}}.$

(H-26)

Since the term on the right-hand side is real, the imaginary part of the left-hand side has to vanish, leading to the final expression

$\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}={\frac {1}{v^{2}(x,y,z)}},$

(H-27)

which is called the eikonal equation. It gives the traveltime value T(x, y, z) for a ray passing through a point (x, y, z) in a medium with velocity v(x, y, z).

A solution to the eikonal equation, T(x, y, z) = constant represents the wavefront at an instant of time. Kinematically, a solution to the eikonal equation (H-27) should also be a solution to the wave equation (H-19). This raises the important question as to under what circumstances the eikonal equation may be considered a good approximation to the wave equation. To investigate this matter, we shall consider a plane wave function as in equation (H-23) but with a spatially varying amplitude P₀(x, y, z) :

$P(x,y,z;t)=P_{0}(x,y,z)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}$

(H-28)

and test its validity as a solution to the wave equation (H-19).

Compute the partial derivative of the wave function P(x, y, z; t) of equation (H-28) with respect to the space variable x

${\frac {\partial P}{\partial x}}=\left({\frac {\partial P_{0}}{\partial x}}+i\omega P_{0}{\frac {\partial T}{\partial x}}\right)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}$

(H-29)

and differentiate the result once more, then simplify to get

${\frac {\partial ^{2}P}{\partial x^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial x^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial x}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+P_{0}{\frac {\partial ^{2}T}{\partial x^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-30a)

Repeat the differentiation for the space variables y

${\frac {\partial ^{2}P}{\partial y^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial y}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+P_{0}{\frac {\partial ^{2}T}{\partial y^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}$

(H-30b)

and z

${\frac {\partial ^{2}P}{\partial z^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial z}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}+P_{0}{\frac {\partial ^{2}T}{\partial z^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-30c)

Next compute the partial derivative of the wave function P(x, y, z; t) of equation (H-28) with respect to the time variable t

${\frac {\partial ^{2}P}{\partial t^{2}}}=-P_{0}\omega ^{2}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.$

(H-30d)

Substitute equations (H-30a,H-30b,H-30c,H-30d) into equation (H-19), and combine the terms into real and imaginary parts to obtain

${\begin{array}{l}\left\{\omega ^{2}P_{0}\left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-\left({\frac {\partial ^{2}P_{0}}{\partial x^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}\right)\right\}\\\quad -i\omega \left[2\left({\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}\right)+P_{0}\left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)\right]={\frac {\omega ^{2}}{v^{2}}}P_{0}.\end{array}}$

(H-31)

Equate the real parts of both sides of equation (H-31) to get

$\left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-{\frac {1}{\omega ^{2}P_{0}}}\left({\frac {\partial ^{2}P_{0}}{\partial x^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}\right)={\frac {1}{v^{2}}},$

(H-32a)

and equate the imaginary parts to get

${\frac {2}{P_{0}}}\left({\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}\right)+\left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)=0.$

(H-32b)

Now we analyze the implications of equations (H-32a,H-32b). To reduce equation (H-32a) to the eikonal equation (H-27), it follows that the second term on the left-hand side of equation (H-32a) must vanish. This is only possible if we let 1/ω go to zero, or in other words, if we make the high-frequency assumption. Thus, for a wave function, such as that given by equation (H-23), with spatially varying amplitudes, the eikonal equation is a good approximation to the wave equation at the high-frequency limit. Additionally, for T(x, y, z) to be a good approximation to the solution of the wave equation, the first term on the left-hand side of equation (H-32b) must be negligible compared to the second term. In fact, the wave function given by equation (H-23) is a special form of the generalized asymptotic ray series solution given by

$P(x,y,z;t)=\sum _{n=0}^{\infty }P_{n}\ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}(i\omega )^{-n},$

(H-33)

where n = 0, 1, 2, …. Note that the terms for n ≥ 1 in equation (H-33) become insignificant for a high-frequency source. Specifically, in the high-frequency limit, the summation given by equation (H-33) reduces to the wave function of equation (H-23).

Since λ = 2πv/ω, where λ is the wavelength, the high-frequency limit is equivalent to small wavelengths. How small the wavelength should be for a solution of the eikonal equation to be a good approximation to the wave equation is an important consideration in the practical validity of the eikonal equation itself. The approximation is valid if the fractional change in the velocity gradient Δv is much less than the frequency v/λ . In practice, this approximation will not be met across a layer boundary with a sharp velocity contrast or in a layer with velocity variations that occur within a spatial extent much less than the wavelength. As a direct consequence of this approximation, the eikonal equation is a good approximation to the wave equation only for a medium in which velocities do not vary rapidly, and especially, do not exhibit sharp discontinuities.

H.3 Finite-difference solution to the eikonal equation

Just as the scalar wave equation itself can be solved using a finite-difference scheme to extrapolate a wavefield from one depth level to another (migration principles), the eikonal equation can be solved using a finite-difference scheme to compute the traveltimes along wavefronts expanding from a source . For simplicity, we shall rewrite the eikonal equation (H-27) for the 2-D case in (x, z) coordinates

$\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}={\frac {1}{v^{2}(x,z)}}.$

(H-34)

Consider the 2-D computation mesh sketched in Figure H-2a. We want to compute the traveltime T of the eikonal equation (H-34) at grid point (x + Δx, z + Δz) using the known traveltimes at grid points (x, z) and (x + Δx, z).

Computing the traveltimes at depth z + Δz from those at depth z means extrapolating T in the z-direction. Rewrite equation (H-34) in the form of an extrapolation equation

${\frac {\partial T}{\partial z}}=\pm {\sqrt {{\frac {1}{v^{2}}}-\left({\frac {\partial T}{\partial x}}\right)^{2}}},$

(H-35a)

which can be rearranged to take the form

$v{\frac {\partial T}{\partial z}}={\sqrt {1-\left(v{\frac {\partial T}{\partial x}}\right)^{2}}}.$

(H-35b)

Note that in equation (H-35b) we also have dropped one of the two solutions and only considered that which yields an increase in traveltime as we march in depth.

Now replace the differential operators with difference operators

${\frac {v(z+\Delta z)T(z+\Delta z)-v(z)T(z)}{\Delta z}}=\left[1-\left({\frac {v(x+\Delta x)T(x+\Delta x)-v(x)T(x)}{\Delta x}}\right)^{2}\right]^{1/2},$

(H-36a)

and rewrite explicitly in terms of T(z + Δz)

$T(z+\Delta z)=T(z)+{\frac {\Delta z}{v(z+\Delta z)}}\left[1-\left({\frac {v(x+\Delta x)T(x+\Delta x)-v(x)T(x)}{\Delta x}}\right)^{2}\right]^{1/2},$

(H-36b)

where we have assumed that v(z)/v(z + Δz) is approximately unity. Also, we have avoided in our algebra crowding the terms in equation (H-36b) with the dual variables (x, z). Now introduce the dual variables explicitly to get the final expression for a finite-difference solution to the 2-D eikonal equation

$T_{i+1}^{k+1}=T_{i+1}^{k}+{\frac {\Delta z}{v_{i+1}^{k+1}}}\left[1-\left({\frac {v_{i+1}^{k}T_{i+1}^{k}-v_{i}^{k}T_{i}^{k}}{\Delta x}}\right)^{2}\right]^{1/2},$

(H-37)

where $T_{i}^{k}$ is the discrete form of T(x, z) as labeled in Figure H-2a.

Figure H-2  Finite-difference mesh used for (a) 2-D and (b) 3-D solution to the eikonal equation. See text for details.

Now consider the 3-D computation mesh sketched in Figure H-2b. We want to compute the traveltime T of the eikonal equation (H-27) at grid point (x + Δx, y, z + Δz) using the known traveltimes at grid points (x, y, z), (x + Δx, y, z) and (x + Δx, y + Δy, z). The 3-D equivalent of equation (H-35b) is

$v{\frac {\partial T}{\partial z}}={\sqrt {1-\left(v{\frac {\partial T}{\partial x}}\right)^{2}-\left(v{\frac {\partial T}{\partial y}}\right)^{2}}}.$

(H-38)

Following the same algebra as for the 2-D case, we obtain the finite-difference solution to the 3-D eikonal equation

${\begin{array}{l}T_{i+1,j}^{k+1}=T_{i+1,j}^{k}+{\frac {\Delta z}{v_{i+1,j}^{k+1}}}\left[1-\left({\frac {v_{i+1,j}^{k}T_{i+1,j}^{k}-v_{i,j}^{k}T_{i,j}^{k}}{\Delta x}}\right)^{2}-\left({\frac {v_{i+1,j+1}^{k}T_{i+1,j+1}^{k}-v_{i+1,j}^{k}T_{i+1,j}^{k}}{\Delta y}}\right)^{2}\right]^{1/2},\end{array}}$

(H-39)

where $T_{i,j}^{k}$ is the discrete form of T(x, y, z) as labeled in Figure H-2b.

Note that in equations (H-36), (H-37), and (H-39), unlike the previous authors , the indices for the velocity v and traveltime T are kept the same. This was done with the stability of the finite-difference solution of the eikonal equation in mind. Actually, Kjartannson has applied the same notion to the finite-difference solution to the scalar wave equation.

See also

• Introduction to earth imaging in depth
• Layer replacement
• 2-D poststack depth migration
• 2-D prestack depth migration
• 3-D poststack depth migration
• 3-D prestack depth migration
• Exercises

External links

find literature about Diffraction and ray theory for wave propagation