Skip to content
BY-NC-ND 3.0 license Open Access Published by De Gruyter Open Access June 27, 2016

Analytical fundamentals of migration in reflection seismics

  • Arnab K. Ray EMAIL logo
From the journal Open Geosciences

Abstract

We consider migration in reflection seismics from a completely analytical perspective. We review the basic geometrical ray-path approach to understanding the subject of migration, and discuss the limitations of this method. We stress the importance of the linear differential wave equation in migration. We also review briefly how a wavefield, travelling with a constant velocity, is extrapolated from the differential wave equation, with the aid of Fourier transforms. Then we present a non-numerical treatment by which we derive an asymptotic solution for both the amplitude and the phase of a planar subsurface wavefield that has a vertical velocity variation. This treatment entails the application of the Wentzel-Kramers-Brillouin approximation, whose self-consistency can be established due to a very slow logarithmic variation of the velocity in the vertical direction, a feature that holds more firmly at increasingly greater subsurface depths. For a planar subsurface wavefield, we also demonstrate an equivalence between two apparently different migration algorithms, namely, the constant-velocity Stolt Migration algorithm and the stationary-phase approximation method.

1 Introduction

In seismic processing, migration of seismic data is a primary requirement. Unless this aspect of processing is addressed satisfactorily, the possibility of making misleading estimates of subsurface depths remains open. Even by a simple quantitative estimate it can be shown that in many cases the error in the spatial measurement of the true reflector position may be of the order of a kilometre [1]. So an accurate measurement of the depth and the local dip of subsurface reflectors is essential. Migration of seismic data, whose principal purpose is to find the correct depths to reflecting interfaces and the corresponding seismic velocities in subsurface rock layers, purports to resolve this difficult question [15]. However, for all its importance, at a basic level migration is usually considered only in passing, under the broad category of reflection seismology, and even so, the understanding of migration is limited to simple geometric ray-path diagrams [3]. Thereafter, in advanced approaches, designed especially for processing geophysicists, migration is taken up with its due degree of mathematical sophistication [1, 2, 4, 5], involving the wave equation, which is a second-order linear partial differential equation, and integral transforms. This, of course, represents a much more comprehensive approach to migration than what could be afforded by a simple geometry-based understanding. Nevertheless, while working under complicated conditions, such as a vertically varying velocity profile of the seismic waves, the results of a wave-field extrapolation are either founded on ad hoc assumptions [1] (on this particular point, the difficulty is acknowledged [2]), or are developed in a manner that appears overwhelming [5]. Further, one often loses sight of the unity of the analytical principles behind various algorithms of migration.

All of these, however, need not be so, as we shall show here. In this paper, we attempt to address certain lacunae in the way seismic migration is understood. The mathematical tools that we implement are quite common in mathematical physics, namely, second-order linear differential equations, Fourier transforms, and the Wentzel-Kramers-Brillouin (henceforth WKB) approximation [68]. To establish how these methods become relevant in migration-related studies, we first discuss the simple geometry-based approach to migration, if only to appreciate its inadequacies (Sections 2 & 3). Then we present a brief review of the use of the wave equation and Fourier transforms [2], which makes the starting point of a more advanced approach (Section 4). Following this, we demonstrate the basic simplicity of the WKB method in deriving a solution for a planar wave with vertical velocity variations (Section 5). And finally, we show that despite apparent dissimilarities in methods, a fundamental unity underlies two separate algorithms on migration (Sections 6 & 7).

Migration is a subject that, in any serious sense, is usually the preserve of a specialist (a processing geophysicist). Bearing in mind the formidable difficulties in this field of study, any adequate treatment of it must have recourse to computer-based techniques [2], and as such is traditionally considered to be beyond a basic introduction to the subject of geophysics [3]. This gap, however, can be bridged through simple but logically consistent analytical methods. We are of the view that the treatment that we have developed in this paper can provide a clear insight into the analytical fundamentals of migration processes. Subsequently, one should feel properly equipped to encounter greater challenges in this complex subject.

2 A basic geometric approach to reflection seismics

In reflection seismics, artificial mechanical waves are generated at predetermined times and spatial positions on the surface of the earth. These waves pass through various earth layers and are reflected from the subsurface interfaces back to the surface of the earth, where they are recorded by receivers. The waves of interest here are compressional (longitudinal) waves, which are also referred to as primary waves [3, 4]. In the simplest possible case, we suppose the surface of the earth to be absolutely flat. On this surface a single source can create waves, that travel through the subsurface, and are subsequently reflected back to the surface from an equally flat single subsurface stratum. This has been schematically represented in the left panel of Fig. 1. We can get the geometrical locus of the physical ray-path in the xt plane, with x defining the source-to-receiver position (the offset), and with the receiver itself gathering reflected signals after an interval of time, t. To this end we apply the Pythagoras theorem to the triangle MDR¯, under the prescription that MR¯=x/2, MD¯=d (the vertical depth of the reflecting interface from the surface of the earth), and DR¯=νt/2 (with v being the constant velocity of the waves propagating through the subsurface). Following some algebraic manipulations, this exercise will deliver the xt locus corresponding to the actual ray-path, known also as the normal moveout equation [3], as

tt02xvt02=1,(1)

which has the form of a hyperbola, with the parameter, t0, defined as t0 = 2d/v. The solution given by Eq. (1) has been geometrically depicted by the continuous curve in the right panel of Fig. 1. From the plot it is clear that different values of d, as a parameter in Eq. (1), will yield a family of hyperbolae (assuming that v retains the same constant value all through). The turning points, the maxima in this case, of this family of hyperbolae, obtained under the condition dt/dx = 0, corresponds to the coordinate (0, t0) in the xt plane. And for any fixed value of x in the xt plane, the entire range of points along the t axis (with t customarily having been scaled along the vertical direction), cutting across the entire family of hyperbolae, will define what is known in seismic processing as a seismic trace [4].

Recast in a slightly different form, Eq. (1) will read as

t=t01+x24d2,(2)

in which t0 can physically be seen to represent the “echo time”, the two-way travel time for the wave, corresponding to a vertical reflection from the subsurface reflector. The term under the square root sign is the normal moveout factor, which appears because the wave reaching the receiver at a distance, x, from the shot point (the location of the source), has not been reflected vertically [3]. Under the condition that the receiver distance, x, is much less than the reflector depth, d, i.e., xd, we can, by a binomial expansion, arrive at an approximate solution for the normal moveout, Δt = t − t0, as

Δtx22ν2t0.(3)

The primary objective here is to determine the depth, d, of the reflector. The vertical echo time, t0, and the normal moveout time, Δt, can be obtained from the reflection data. The corresponding receiver distance, x, can also be known from the data. Taken together, all of these allow us to determine the value of v from Eq. (3). This, along with the definition of t0, will also enable us to determine d. Obviously, a precise extraction of the value of v from the reflection data is necessary for determining the reflector depth accurately. Knowing the velocity correctly is equally important when we consider more complicated instances of reflection seismology, involving non-horizontal reflectors of continually varying gradients. Besides this, with v being dependent on the elastic properties of the material through which the wave propagates [4, 9], getting a correct measure of the velocity also informs us of the true nature of the subsurface material.

Figure 1 A geometric impression pertaining to a single source-receiver pair and a single horizontal reflector(a) A schematic representation of reflection seismics from a flat subsurface stratum, using a single source-receiver pair. The waves originate from the source position, S. They are reflected from an absolutely flat interface in the subsurface, back to the receiver position, R. The midpoint between S and R, shown as M, is vertically above the point from where the reflection occurs, indicated as D. For many source-receiver arrays placed as symmetrically about M, as it has been shown in the figure, the reflection point D is a common depth point for all such source-receiver pairs. The source-to-receiver distance, SR¯$\overline {{\rm{SR}}} $, is x. The depth to the reflecting interface, MD¯$\overline {{\rm{MD}}} $, is given as d. The signal takes a time, t, travelling with a constant velocity, v, to reach R from S, after undergoing a reflection at D. Triangle SDR¯$\overline {{\rm{SDR}}} $ is an isosceles triangle, while triangle MDR¯$\overline {{\rm{MDR}}} $ is a right-angled triangle.(b) The locus of the ray-path between a source and a receiver in the x–t plane. In this plot, t increases in the downward direction along the vertical axis. The origin of coordinates, (0, 0), is at O. The continuous curve corresponds to the source-receiver array shown in the left panel of this figure. The peak of this curve is at the point P, whose coordinates are (0, t0) in this x–t plot. The dotted curve corresponds to the situation shown in the left panel of Fig. 2. The two horizontal lines are the tangents passing through the maxima of the two hyperbolae. It is evident from this plot that inclined subsurface reflectors cause a shift in the locus of the ray-path. In this case the maximum point on the curve has shifted to the point P1. With continuous change in the dip angle of the subsurface reflector, the path traced out by the peaks of the hyperbolae is an ellipse, expressed mathematically by Eq. (5).
Figure 1

A geometric impression pertaining to a single source-receiver pair and a single horizontal reflector

(a) A schematic representation of reflection seismics from a flat subsurface stratum, using a single source-receiver pair. The waves originate from the source position, S. They are reflected from an absolutely flat interface in the subsurface, back to the receiver position, R. The midpoint between S and R, shown as M, is vertically above the point from where the reflection occurs, indicated as D. For many source-receiver arrays placed as symmetrically about M, as it has been shown in the figure, the reflection point D is a common depth point for all such source-receiver pairs. The source-to-receiver distance, SR¯, is x. The depth to the reflecting interface, MD¯, is given as d. The signal takes a time, t, travelling with a constant velocity, v, to reach R from S, after undergoing a reflection at D. Triangle SDR¯ is an isosceles triangle, while triangle MDR¯ is a right-angled triangle.

(b) The locus of the ray-path between a source and a receiver in the xt plane. In this plot, t increases in the downward direction along the vertical axis. The origin of coordinates, (0, 0), is at O. The continuous curve corresponds to the source-receiver array shown in the left panel of this figure. The peak of this curve is at the point P, whose coordinates are (0, t0) in this xt plot. The dotted curve corresponds to the situation shown in the left panel of Fig. 2. The two horizontal lines are the tangents passing through the maxima of the two hyperbolae. It is evident from this plot that inclined subsurface reflectors cause a shift in the locus of the ray-path. In this case the maximum point on the curve has shifted to the point P1. With continuous change in the dip angle of the subsurface reflector, the path traced out by the peaks of the hyperbolae is an ellipse, expressed mathematically by Eq. (5).

3 Dipping reflectors and a case for migration

It is scarcely true that the layering of actual geological strata will conform to the neat and orderly picture implied by the diagram in the left panel of Fig. 1. If anything, not only will the subsurface strata be quite great in number, but also each stratum will have locally varying spatial gradients. This will cause much difficulty in devising an accurate theoretical model of the physics of reflection seismics, a fact that can be illustrated even by considering the very simple case of a single dipping interface with a constant spatial gradient.

When the reflecting interface is inclined at a constant angle, θ, with respect to the horizontal, as it has been shown in the left panel of Fig. 2, we have to make suitable modifications in the analysis presented in Section 2, to retain the canonical form of the hyperbolic normal moveout equation. This can be achieved by noting that the perpendicular distance from the source to the inclined reflector is SN¯=d, the source-to-receiver distance is SR¯=x, and the total travel path is SD1¯+D1R¯=νt. With this information, even upon accounting for the dip in the reflector, the invariant canonical form of the hyperbola can be obtained as

tt¯02x¯νt¯02=1,(4)

with x¯ and t¯0 being defined by the transformations, x¯=x+2dsinθ and t¯0=(2dcosθ)/ν, respectively. From Eq. (4), the coordinates of the maximum point of the hyperbola (xm, tm) in the xt plane (with t increasing in the downward direction) will be given by the condition dt/dx = 0. This will place the peak of the hyperbola at (2dsinθ,t¯0). Going back to the case of the horizontal reflector (θ = 0), given by Eq. (1), the peak of the hyperbola here is to be seen at (0, t0). Therefore, a distinct shift is seen to arise because of the dip of the subsurface reflector, and this shift has been shown by the position of the dotted hyperbola in the right panel of Fig. 1. The locus of this shift can also be traced on the xt plane. Noting that the coordinates of the maximum point of the normal moveout hyperbola are given by Xm = −2d sin θ and tm=t¯0, to, we could, by making use of the standard trigonometric result, sin2θ + cos2θ = 1, derive the condition,

xm2d2+tm2d/ν2=1.(5)
Figure 2 A geometric view of the complications due to an inclined subsurface reflector(a) Reflection from an inclined subsurface reflector, using a single source-receiver pair placed symmetrically about the midpoint, M. The reflector is at a dip angle, θ, with respect to the horizontal. The actual point of subsurface reflection, D1, is not vertically below M (unlike the point D), and its location varies for each set of source-receiver pair placed symmetrically about M. There is no common depth point for these multiple source-receiver pairs. The perpendicular distance from the source to the reflector is SN¯=d$\overline {{\rm{SN}}}  = d$. The total travel path is SD1¯+D1R¯=νt$\overline {{\rm{S}}{{\rm{D}}_{\rm{1}}}}  + \overline {{{\rm{D}}_{\rm{1}}}{\rm{R}}}  = \nu t$, and the source-to-receiver distance, SR¯=x$\overline {{\rm{SR}}}  = x$.(b) For an inclined reflector, imaging the depth vertically from the surface of the earth, gives a false position of the reflector. The actual point of reflection on the inclined plane is indicated by N, but a vertical downward projection from S will give a wrong impression that the reflection point is at D. The line, SN¯$\overline {{\rm{SN}}} $, is normal to the reflector interface. Migration gives the correct position of the true reflector interface. From a practical point of view, the inclined reflector never has a simple and fixed linear profile, but instead has local spatial variations. This contributes further errors in the measurements.
Figure 2

A geometric view of the complications due to an inclined subsurface reflector

(a) Reflection from an inclined subsurface reflector, using a single source-receiver pair placed symmetrically about the midpoint, M. The reflector is at a dip angle, θ, with respect to the horizontal. The actual point of subsurface reflection, D1, is not vertically below M (unlike the point D), and its location varies for each set of source-receiver pair placed symmetrically about M. There is no common depth point for these multiple source-receiver pairs. The perpendicular distance from the source to the reflector is SN¯=d. The total travel path is SD1¯+D1R¯=νt, and the source-to-receiver distance, SR¯=x.

(b) For an inclined reflector, imaging the depth vertically from the surface of the earth, gives a false position of the reflector. The actual point of reflection on the inclined plane is indicated by N, but a vertical downward projection from S will give a wrong impression that the reflection point is at D. The line, SN¯, is normal to the reflector interface. Migration gives the correct position of the true reflector interface. From a practical point of view, the inclined reflector never has a simple and fixed linear profile, but instead has local spatial variations. This contributes further errors in the measurements.

The form of this conic section is that of an ellipse. It can now be concluded that the peak of the normal moveout curve will trace an elliptical trajectory in the xt plane, with a continuous change in the dip angle, θ.

While this entire algebraic exercise looks very elegant for pedagogic purposes, it fails hopelessly in addressing a serious practical problem. The derivation of Eq. (4) has been done with the help of a single source-receiver pair. It assumes that every source-receiver pair in a full array of sources and receivers has a common point of reflection from the inclined subsurface. A priori this is not possible to achieve in practice. If anything, for a realistic arrangement of source-receiver arrays, each source-receiver pair has its own particular reflection point. Therefore, unlike the simple representation in the left panel of Fig. 1, there is no common depth point vertically below the common midpoint of the source-receiver arrays. Overlooking this crucial point will convey an erroneous impression of the subsurface depth. So this fact precludes our simple approach of determining the value of d, with the arguments that we have presented so far. Our difficulties, however, do not end here. In actual fact the situation is vastly more complicated. The dipping reflectors have a varying gradient from one spatial position to another, as opposed to being at a constant incline with respect to the horizontal. Besides this, there are geological faults, which in general are threedimensional inhomogeneities of various physical properties of the earth’s crust, such as its density and electrical conductivity [10]. These physical properties undergo variation at the faults [10, 11], just as the propagation velocity of seismic waves also has a variation with respect to the depth. Taken together, all of these facts render the accurate imaging of the subsurface quite a formidable task. Nevertheless, all reflection seismic records must be corrected for non-vertical reflections. This evidently complicated process of correction is what is implied by the term migration. It purports to shift all dipping reflections to their true subsurface positions [1, 3]. We need to stress another important point here. Till this point, the entire discussion has accounted for dip and subsurface variations in the inline direction [1], i.e., the direction along which the source-receiver arrays have been lined up. This means that the analytical treatment presented so far has been confined to a 2D plane only — along the vertical depth (given by the z coordinate) and along the offset direction (given by the x coordinate). However, variations in the subsurface will practically also occur in the crossline direction [1], i.e., along the y axis. So a comprehensive migration process must account for variations along all the three spatial coordinates, i.e., it should be a complete 3D process, accounting for changes in the gradient of a subsurface reflector.

4 Migration with the differential wave equation

By now we have fully appreciated the limitations of a simple geometric approach to migration, and, along with it we have also argued for the importance of a 3D migration process. To carry out such a migration effectively, we have to study the propagation of a compressional wavefield, P = P(x, y, z, t). This wavefield is described locally by the linear differential wave equation,

2P1ν22Pt2=0,(6)

in which ▽2 is the standard Laplacian operator [6], and in which v has a dependence on all the three spatial coordinates. Expressed in terms of individual spatial components, this second-order partial differential equation in both space and time is

2x2+2y2+2z21ν22t2P=0.(7)

On the surface of the earth the upcoming seismic wavefield is measured as P(x, y, 0, t). For fixed values of x and y, the signal is a function of t only, giving a single trace. Extended over the whole xy plane, this will build an entire 3D volume of seismic traces. Given this boundary condition on the surface of the earth, we have to find the wavefield, P(x, y, z, 0), which is actually the reflectivity of the earth in its recesses, and which also defines an initial condition at t = 0. Thus, this process requires projecting the surface wavefield to a depth z, corresponding to t = 0. This is the mathematical essence of migration [1, 2].

Solving Eq. (7) is no easy feat, because of its dependence on four independent variables. However, under the simple working assumption that v is a constant, it is possible to go a long way, by adopting the method of Fourier transforms [6]. By definition, the Fourier transform of the wavefield P(x, y, z, t) over the coordinates x, y and t is given as [6]

P~(kx,ky,z,ω)=1(2π)3/2P(x,y,z,t)exp(ikxx+ikyyiωt)dxdydt.(8)

In terms of this transformed wavefield, we can recast Eq. (7) as a second-order ordinary differential equation, which is

d2P~dz2+ω2ν2kx2ky2P~=0,(9)

with P~ now being understood to be P~=¯P~(kx,ky,z,ω). It is worth stressing the importance of the method that has led to Eq. (9). The starting point was Eq. (6), which is a second-order partial differential equation in four variables. The Fourier transform on x, y and t into their respective conjugate variables, kx, ky and ω, has converted this partial differential equation to a simple second-order ordinary differential equation. For mathematical ease, this is the most important purpose that a Fourier transform can serve. Of course, working in the Fourier domain has other practical advantages, especially from the perspective of the processing geophysicist, like the elimination of spikes and noise from seismic data [4].

Following the Fourier transform, Eq. (9) can now be applied to extrapolate the wavefield from its value on the surface of the earth (at z = 0) to the required depth (where z ≠ 0). Working under the simplifying assumption that v is constant [1], it becomes quite easy to solve Eq. (9), which, in fact, has the form of the simple harmonic oscillator [6]. Imposing the definition [1, 2]

kz2=ω2ν2kx2ky2,(10)

which is a dispersion relatio[1], we can proceed to derive an integral solution of Eq. (9) as

P~(kx,ky,z,ω)=P~(kx,ky,0,ω)exp(ikzz).(11)

Actually, the most general form of Eq. (11) is given by the superposition of an upcoming wave and a downgoing wave, i.e., exp(±ikzz). However, only the solution with the negative sign need be retained here, since it corresponds physically to upcoming waves, which, for our present purpose, are the relevant modes [1]. This whole analysis shows how it should be possible to extrapolate the surface wavefield to a given depth below the surface of the earth. While it is quite a simple task to arrive at this result under the assumption of constant velocity, we can also derive an asymptotic solution for the extrapolation exercise when the velocity is dependent on the depth, i.e., v = v(z). This treatment has been presented in Section 5, but at this stage the one primary fact that stands out is that so far as migration is concerned, the result given by Eq. (11), i.e., extrapolating the surface data to the depth, represents a crucial step. Beyond this point, there can be many variations in the subsequent execution of the migration process, depending on varying contingencies. In what follows, we shall dwell upon the principles of some of these methods, and derive meaningful conclusions from them.

5 Extrapolating the wavefield with depth-dependent velocity: The WKB approximation

We have seen how easily Eq. (9) can be integrated when v is constant. Supposing now that v has a dependence on the z coordinate, i.e. v ≡ v(z), it will still be possible to derive an asymptotic solution of Eq. (9), with the help of the WKB approximation [7, 8]. We offer this treatment as a simple and a more elucidating alternative to the way this mathematical problem is treated even in very advanced texts on migration [1]. We work with the assumption of a planar wavefront propagating downwards with a velocity, v(z). Conditions pertaining to this restriction on the propagating wave would be achieved satisfactorily, if there should be a large number of sources placed on the (x, y, 0) plane. In this situation, the effective wavefront propagating downwards will be planar, except at the fringes. The error due to the fringe effect can be reduced by setting up a very wide field of the seismic survey[2].

Now, since v has a dependence on z, so should kz have a dependence on z, by virtue of Eq. (10). Therefore, we can write

d2P~dz2+kz2(z)P~=0,(12)

for which, a trial solution in the most general form can be set down as

P~=A(z)eiα(z),(13)

with A being an amplitude function and α being a phase function. Substituting this trial solution in Eq. (12) gives

1Ad2Adz2dαdz2+kz2(z)+idαdzddzIndαdzA2=0.(14)

The real and imaginary components above are set to zero separately. Thereafter, a solution of the imaginary part is extracted, and it reads as

A=Cdαdz1/2,(15)

where C is a constant of integration. The foregoing result indicates that knowing either A or α independently suffices to give the full solution implied by Eq. (13).

So far the mathematical treatment has been exact. To derive a solution from the real component of Eq. (14), it is necessary now to apply the WKB method. We note that in studies on migration (see [5] and references therein), this method is regularly implemented to account for the vertical variation of the velocity function, but often such results appear inaccessible. Therefore, our purpose here is to establish and appreciate the simple analytical essence of the WKB method, which is commonly invoked in fundamental physical problems dealing with the propagation of waves. To this end we approximate the amplitude function, A(z), in Eq. (13), to vary much more slowly over z than the phase function, α(z). Consequently we can neglect the second derivative of A with respect to z in the real part of Eq. (14). We point out that this approximation is much more selfconsistent than completely ignoring the spatial gradient of the velocity function, as it is occasionally prescribed [1]. So with our chosen approximation, we are led to the phase solution,

α(z)kz(z)dz,(16)

(taking the more practically relevant negative solution, as usual) which, used along with Eq. (15), gives an aysmptotic solution that reads as

P~Ckzexpikz(z)dz.(17)

The simplicity and the compactness of the foregoing result are now plainly visible. Through a series of logically consistent steps, we have derived an asymptotic solution for the downwardly propagating wavefield, as given by Eq. (17), which bears analytical expressions of both the amplitude and the phase. We may ask how the use of the WKB approximation can be justified self-consistently. On using the solutions given by Eqs. (15) and (16), it can be seen that the second derivative of A is given as

d2Adz2A2ωk2ν2252ωk2ν2×ddz(Inν)2d2dz2(Inν).(18)

On the right hand side of Eq. (18) there are derivatives of the logarithm of v. Now, the logarithm of a function varies much more slowly than the function itself. This reasoning is applied to ln v in comparison to v itself. Therefore, we can argue that with the second derivative of A depending on the derivatives of ln v, the requirement of A being a slowly varying function of z is satisfactorily fulfilled, and we may safely neglect the term containing the second derivative of A in Eq. (14). Indeed, this is all the more appropriate at greater depths, where, with v being an increasing function of z, the logarithmic terms strongly temper the variation of A arising due to the velocity variations. So the asymptotic solution given by Eq. (17) has increased validity at greater depths.

6 The 3D phase-shift migration in the frequency-wavenumber domain

We discussed in Section 4 that the definition of migration necessitates projecting seismic data from the boundary condition, z = 0 for all t, to the initial condition, t = 0 for all z. To this end, we find it expedient to look first at the inverse Fourier transform corresponding to Eq. (8). This is

P(x,y,z,t)=1(2π)3/2P~(kx,ky,z,ω)exp(ikxxikyy+iωt)dkxdkydω.(19)

We also know by now how it is possible to relate P(kx, ky, z, ω) to the surface seismic data with the help of Eq. (11). Using this condition in Eq. (19) and setting t = 0, we get

P(x,y,z,0)=1(2π)3/2P~(kx,ky,0,ω)exp(ikxxikyy+ikzz)dkxdkydω,(20)

which is the equation of the phase-shift method, suited for 3D migration [1]. We can substitute ω in this equation, by making use of Eq. (10), under the constraint that v is constant, and kx and ky are unchanged [1]. This delivers the relation

dω=νkzkx2+ky2+kz2dkz,(21)

leading ultimately to

P(x,y,z,0)=1(2π)3/2νkzkx2+ky2+kz2P~(kx,ky,0,νkx2+ky2+kz2)exp(ikxxikyy+ikzz)dkxdkydkz,(22)

which is the equation of the constant-velocity 3D Stolt migration [1]. The derivation of the result in Eq. (22) is certainly consistent with the mathematical understanding of migration in Section 4. However, short of resorting to an involved numerical treatment, the full 3D problem does not lend itself very easily to deriving any analytical insight. Hence, for simplicity of understanding, we consider a very special case of this migration method, for a wavefield that has a spatial dependence on the z coordinate only, which in effect becomes a vertically propagating planar wavefield. For such a wavefield, with no dependence on the x and y coordinates, and propagating in a planar front along the z axis only, we can set down the wavefield as P = P(z, t). In terms of the Fourier inverse transform it can be expressed as

P(z,t)=1(2π)1/2P~(z,ω)exp(iωt)dω.(23)

Now, we know already that Eq. (11) hints at an extrapolation relation going as P~(z,ω)=P~(0,ω)expikzz, for the special case we are considering here. From Eq. (10) we also get kz = ω/v, since there is no involvement of the x and y coordinates. Use of all these conditions in Eq. (23) results in

P(z,t)=1(2π)1/2P~(0,ω)expiωzν+tdω,(24)

which, for t = 0, gives

P(z,t)=1(2π)1/2P~(0,ω)expiωzνdω=P(0,z/ν).(25)

We also know that t and ω are conjugate variables in a Fourier transform. Their mathematical connection is given by

P(0,t)=1(2π)1/2P~(0,ω)expiωdω,(26)

which can then be compared with what is given in Eq. (25). This immediately yields

z=νt,(27)

which is a result that connects the depth of a reflecting stratum with the time of the reception of a signal reflected from that stratum (under the simple condition that the velocity of the wave propagation is unchanging). An important point to note here is that the simple appearance of Eq. (27) has been possible because the plane wave accommodates only a vertically propagating mode in Eq. (10), and so by default this mode has been decoupled from all the other modes. This decoupling should be impossible to achieve for higher spatial dimensions, as a look at Eq. (10) ought to reveal.

7 The stationary-phase approximation

We now derive the traveltime equation for the 3D zero-offset case, using the method of the stationary-phase approximation. From Eq. (19) we can, after using the extrapolation condition given by Eq. (11), write a compact relation that reads as

P(x,y,z,t)=1(2π)3/2P~(kx,ky,0,ω)exp(iϕz)dkxdkydω,(28)

in which

ϕ=kzkxxzkyyz+ωtz.(29)

It is easy to see that the integrand in Eq. (28) has an amplitude part and a phase part. If the phase is to vary much more rapidly than the amplitude (a requirement that is consistent with the WKB analysis presented in Section 5), then the rapid oscillations over most of the range of the integration will result in an average value of nearly zero. The exception to this principle occurs when Φ is stationary [6, 7]. This (the stationary-phase approximation) is quantified by δΦ = 0, which is a condition for the extremum value of Φ. The major contribution to the integral, therefore, comes from the extremum points of the phase function, Φ. This condition helps in identifying the physical coordinates from where we obtain the most significant contribution to the signal received on the surface of the earth.

Going by Eq. (10), we know that kz has a dependence on kx, ky and ω. Taking this in conjunction with Eq. (29), and noting that the relevant variables for the integration in Eq. (28) are kx, ky and ω, the required extremum condition for Φ can consequently be stated as

δϕ=ϕkxkx+ϕkyky+ϕωω=0.(30)

The result above can only mean that all the three partial derivatives in Eq. (30) are to be set individually to zero. This will give us three distinct mathematical conditions. Using Eq. (29) we first get

ϕkx=kzkxxz=0,(31)

after which we go back to Eq. (10) and obtain

kzkx=kxkz,(32)

a result that we then combine with Eq. (31) to finally get

zkx=kzx.(33)

Working similarly with ky in Eqs. (29) and (10), we get

zky=kzy,(34)

and likewise with ω, we get

zων=νkzt.(35)

Now we square Eqs. (33), (34) and (35). Then we subtract the squares of Eqs. (33) and (34) from the square of Eq. (35). Thereafter, making use of Eq. (10) once again, we derive the traveltime equation for the 3D zero-offset case as

x2+y2+z2=(νt)2.(36)

This defines a sphere of radius vt in the x–y–z coordinate space, and a hyperboloid in the x–y–t space for a constant value of z. Considering only the inline direction and the depth (i.e., when y = 0), Eq. (36) will define a semicircle below the surface of the earth in the x–z plane, and the usual hyperbola in the x–t plane. However, a much more interesting result is observed for the simplest possible case of the planar wavefront, for which the spatial dependence is only on the z coordinate. In that case, the result in Eq. (36) will converge simply to the linear solution z = vt, which is exactly what Eq. (27) also furnishes. So effectively this simple special case establishes a mathematical equivalence between the stationary-phase approximation method and the constant-velocity Stolt migration algorithm.

8 Summary and concluding remarks

This study is on the mathematical fundamentals of migration in reflection seismics. The whole discussion is to be viewed under two broad categories. The first deals with a basic review of migration and its associated difficulties. The second discusses some analytical methods applied to the field of migration with some consequent insights. The main outcomes of this study are:

  1. The implementation of the WKB approximation on a planar wave with a vertical velocity variation, resulting in a compact analytical solution for both the amplitude and the phase of the wave.

  2. The equivalence between the constant-velocity Stolt migration algorithm, and the stationary-phase approximation for a downwardly propagating planar wave.

Our study is designed to contribute to the fundamental understanding of the subject of migration. However, at the end we must mention that beyond a basic understanding of the subject, the actual implementation of a migration algorithm, in any realistic sense, is a formidable task.

None of the migration algorithms implemented in the industry today gives a complete prescription for handling all dips and velocity variations. Looking at some of the algorithms we have discussed here, the ones based on the integral solution of the wave equation can become cumbersome when they encounter lateral velocity variations [1]. And frequency-wavenumber algorithms are effective against lateral velocity variations only to a limited extent [1].

Acknowledgements

The author thanks J. K. Bhattacharjee, B. Bleines, S. Roy Chowdhury and S. Sarkar for their comments. Thanks are also due to A. R. Dhakulkar, A. Kumar, T. Mukherjee and P. Pathak for providing help in some technical matters. A careful reading of the manuscript by P. J. Wiita is gratefully acknowledged.

References

[1] Yilmaz Ö., Seismic Data Analysis. Society of Exploration Geophysicists, Tulsa, Oklahoma, 200110.1190/1.9781560801580Search in Google Scholar

[2] Claerbout J. F., Imaging the Earth’s Interior. Stanford University, California, 1985Search in Google Scholar

[3] Lowrie W., Fundamentals of Geophysics. Cambridge University Press, Cambridge, 1997Search in Google Scholar

[4] Mari J.-L., Glangeaud F., Coppens F., Signal Processing for Geologists and Geophysicists. Éditions Technip, Paris, 199910.2516/ifpen/2011002Search in Google Scholar

[5] Stolt R. H., Weglein A. B., Seismic Imaging and Inversion: Application to Linear Inverse Theory. Cambridge University Press, Cambridge, 201210.2172/1052406Search in Google Scholar

[6] Arfken G. B., Weber H. J., Mathematical Methods for Physicists. Academic Press, San Diego, 2001Search in Google Scholar

[7] Mathews J., Walker R. L., Mathematical Methods of Physics. Pearson Education, Singapore, 1970Search in Google Scholar

[8] King A. C., Bilingham J., Otto S. R., Differential Equations. Cambridge University Press, Cambridge, 200310.1017/CBO9780511755293Search in Google Scholar

[9] Landau L. D., Lifshitz E. M., Theory of Elasticity. Butterworth-Heinemann, Oxford, 1986Search in Google Scholar

[10] Sarlis N., Lazaridou M., Kapiris P., Varotsos P., Numerical model of the selectivity effect and the ΔV/L criterion. Geophysical Research Letters, 1999, 26, 3245-324810.1029/1998GL005265Search in Google Scholar

[11] Varotsos P., Sarlis N., Lazaridou M., Kapiris P., Transmission of stress induced electric signals in dielectric media. Journal of Applied Physics, 1998, 83, 60-7010.1063/1.366702Search in Google Scholar

Received: 2015-6-17
Accepted: 2015-10-30
Published Online: 2016-6-27
Published in Print: 2016-6-1

© 2016 A. K. Ray, published by De Gruyter Open

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.

Downloaded on 2.3.2024 from https://www.degruyter.com/document/doi/10.1515/geo-2016-0025/html
Scroll to top button