Part II: Attaining accuracy with depth-variable velocity dip moveout

Nowadays, a standard industrial seismic processing almost always involves the constant-velocity dip moveout correction. Two main features make the process attractive. First, the computational cost is low compared to prestack migration in both two-dimensional and three-dimensional processings. Secondly, because constant-velocity dip moveout is ``independent of velocity'', it may come before velocity analysis, removing the effects of dip ().

However, the hypothesis of constant velocity is somewhat inconsistent with the concept of velocity analysis which results in the construction of a depth-variable velocity model. Recently, several methods for depth-variable velocity dip moveout (, , ) undeniably improved the zero-offset stack section. However, one can wonder if they clearly improve the result of the post-stack migrated section.

How variable velocity dip moveout improves post-stack migration

In this section I show that the depth-variable velocity dip moveout strongly differs from the constant-velocity dip moveout not by the shape of the operator (although the operator is squeezed horizontally when the velocity increases with depth (, )) but by the amplitude distribution along the operator. In the second section, a synthetic data example shows how the processing flow NMO-v(z) DMO-post-stack migration is comparable to pre-stack migration.

Some amplitude issues

Jakubowicz gave an elegant formulation of f-k dip moveout where a finite number of dips are processed separately and stacked according to equation (11),  

$$P_{\rm DMO} = \sum_i \Delta D_i \ P_{D_i},$$ (11)

where P_Di is the wavefield, NMO-corrected for a single dip D_i, and $\Delta D_i$ is the width of the dip-range surrounding D_i. Obviously, if the dip sampling is irregular, the dip-decomposed wavefields stack with a different weight. Beasley derived an amplitude scheme for Kirchhoff dip moveout from Jakubowicz's weighting by using the following practical approximation:  

$$\Delta D_i \approx \frac{\partial^2 t_0}{\partial x^2} \Delta x.$$ (12)

To the first order, the amplitude along the dip moveout operator is proportional to its time curvature. This approximation is close to the amplitude scheme mathematically derived by Black et al. from (f,k) dip moveout (). Thus, Jacubowicz's dip moveout by dip-decomposition, Black's Kirchhoff dip moveout, and (f,k) dip moveout have a very similar amplitude distribution.

In the general case of depth-variable velocity, there is no analytic expression of the amplitude along the operator but the dip-decomposition method gives us a qualitative idea. Indeed, the operator of variable-velocity dip moveout may be build point by point for a range of dips regularly sampled. Then, a region where the points are close to each other allows a significant stack of the dipping segments of the operator, and thus corresponds to an area of high amplitude. Figure (imprs) shows an impulse response of dip moveout in a time-variable velocity model whose profile is represented in Figure (vel).

 

vel Figure 16 Interval velocity model.

 

imprs
Figure 17 Impulse response of v(z) 2-D DMO for an offset of .8 km and a total two-way traveltime of 1. s. The stars along the operator correspond to a uniform sampling of dips. The denser are the points, the higher is the amplitude of the operator. The dashed line represents the elliptic support of the constant-velocity DMO impulse response for the same offset and traveltime. Both the solid and dashed lines have a similar shape in the region of high amplitudes, differing only by the triplications along the solid line. The unrealistic points beyond the offset show that there is a problem of convergence of the v(z) DMO algorithm for steep dips.

Artley's v(z) dip moveout method which produced Figure (imprs), requires a non-linear inversion of a $4 \times 4$ system of equations for each time, offset, and dip we consider. Though it is a time consuming process, it provides us with amplitude information. First, the triplication of the impulse response have a low energy for this velocity model (a trough in the velocity profile would show a high amplitude triplication). Secondly, the amplitude varies a lot along the operator, starting high at gentle dips, decreasing, and increasing again just before a triplication.

The shape of a v(z) DMO operator differs from the constant-velocity DMO operator essentially by a number of triplications along its branches (Figure (imprs)). However, since these triplications are generally low-amplitude, the shape is not the most discriminative feature of the operator. The amplitude variations along the operator strongly depend on the velocity model, and therefore play an important role in the dip moveout correction.

Comparison with constant-velocity dip moveout

There is no doubt that replacing the constant-velocity dip moveout step by a variable-velocity dip moveout step improves the zero-offset section after stack. However, the final result of a standard processing flow is the post-stack migrated section. Thus, the idea of comparing the pre-stack migrated section to the post-stack migrated sections with and without the constant velocity assumption is natural. David Lumley provided the synthetic data example that appears in the next three figures. This synthetic model includes four diffractors at various depth and midpoint positions, two horizontal events, and a 45-degree dipping event, overlaid by a strong, constant RMS-velocity gradient ($V_{\rm RMS}^{({\rm km/s})} = 1.5 + 2. \tau^{({\rm s})}$).

Results of constant-velocity DMO

The processing flow involving a constant-velocity dip moveout step is displayed in Figure (cvmig). The three deepest diffractors are unfocused whereas the shallow diffractor is focused but low-amplitude. The horizontal reflectors are well imaged but the dipping event has a lower amplitude than expected though it is not mispositioned.

 

cvmig
Figure 18 Post-stack migrated section after NMO, constant-velocity DMO and stack. The deepest diffractors are unfocused and the dipping event is low-amplitude.

Results of approximate variable-velocity DMO

Castle proposes a method conceptually similar to Hale's squeezed dip moveout operator but different in its implementation. Instead of applying a standard normal moveout correction ($t_n^2 = t^2 - x^2 / v_{\rm RMS}^2$)followed by a dip moveout correction with a squeezed operator, the method uses de Bazelaire's normal moveout correction with shifted hyperbolae:

$$t_p^2 = (t-t_n+t_p)^2 - \frac{x^2}{\nu^2},$$ (13)

where t_p = t_n/S, $\nu^2 = S \mu_2$, S is the squeezing factor given by $S = \mu_4 / \mu_2^2$, and $\mu_k$ is the order-k momentum ($\mu_2 = v_{\rm RMS}^2$). Castle shows that applying the constant-velocity dip moveout operator in the t_p time domain and then shifting back the data according to the equation

t₀ = t_p S

is equivalent to squeezing the dip moveout operator. Both methods approximate the exact dip moveout operator for depth-variable velocity by wiping out the triplication apparent in Figure (imprs).

Figure (svmig) results from post-stack migration after the approximate dip moveout step. The four diffractor points are better focused with this flow than with the flow that involves constant-velocity dip moveout. Unfortunately, the method, apart from approximating the shape, also simplifies the amplitude distribution along the operator. The dull appearance of the reflectors is a consequence of this improper handling of the amplitudes.

 

svmig
Figure 19 Post-stack migrated section after NMO, approximate variable-velocity DMO and stack. The diffractors and the dipping event are better focused than in Figure (cvmig) but have unbalanced amplitudes.

Result of variable-velocity DMO

For this processing flow, I use the variable-velocity DMO algorithm developed by Artley . The method, accurate for any depth-variable velocity media, solves a system of equation accounting for the location and the dip of any reflector point. An extension of this system to three dimensions is proposed in the next section (page ). The similarity between the pre-stack migration (Figure (psmig)) and the post-stack migration after depth-variable velocity dip moveout Figure (vvmig) is not surprising because the two processing are equivalent when there is no lateral velocity variation. However, the clear improvement with respect to the standard processing makes the variable-velocity dip moveout a very useful tool.

 

vvmig
Figure 20 Post-stack migrated section after NMO, exact variable-velocity DMO and stack. The diffractors and the dipping event are well imaged.

 

psmig
Figure 21 Pre-stack migration. Click on the following button to see a movie of the four previous pictures.

Conclusion

The main difference between constant-velocity and variable-velocity dip moveout is not the shape of the operator but the amplitude distribution. Because this distribution highly depends on the velocity model, variable-velocity dip moveout considerably improves the post-stack migration. However, the drawback is its high computational cost. For example, Artley's method is more costly than a pre-stack migration, because the algorithm does much more than compute the zero-offset traveltime: it performs a ray-based pre-stack migration and uses the computed zero-offset traveltime to apply the dip moveout correction. Because both Meinardus's and Artley's algorithms are based on the dip-decomposition idea, it is possible to reduce the computational cost by cutting down the number of dips to process. On the other hand, Castle's method differs from the constant-velocity dip moveout flow by involving an extra time-variable shift of the data. At a small additional computational cost, the approximate dip moveout flow improves the quality of imaging after post-stack migration. Furthermore, this method can be adapted in three dimensions by using the constant-velocity algorithm described in the first part (page ).

Apart from all speed considerations, a precise dip moveout method requires both assumptions of a depth-variable velocity and a three-dimensional model. The next section explains how Artley's two-dimensional depth-variable velocity dip moveout method can be extended to three dimensions.

Equations for three-dimensional dip moveout
in depth-variable velocity media

  Artley introduced an original method to perform v(z) dip moveout in a two-dimensional earth model. The process uses ray tracing tables and solves a system of equations accounting for the location and the dip of the reflection point. This section describes an extension of Artley's method and derives a new set of equations for the 3-D case.

The ray parameter in three dimensions

In a two-dimensional earth model, the ray parameter is given at any point of the ray by the relationship  

$$p = \frac{\sin\theta}{v} ,$$ (15)

where $\theta$ is the inclination of the ray path with respect to the vertical axis and v is the local velocity. We can see that the ray parameter is simply the projection length of the ray path vector $\bf r$($\Vert{\bf r}\Vert = r = 1/v$) on the earth surface. In a 3-D v(z) model of the earth, the rays travel in a vertical plane. In this case, the ray parameter, being the projection of the ray path vector, is a two-dimensional vector that can be expressed in either cartesian or polar coordinates, as follows:  

$${\bf p} = \left( \begin{array} {c} p_x \\ p_y \end{array}... ...\begin{array} {c} \cos\phi \\ \sin\phi \end{array} \right) ,$$ (16)

where $\phi$ is the strike of the vertical plane containing the ray and $p = \sqrt{ p_x^2 + p_y^2 }$.

Deriving the system of equations

Vectorial relationships on the earth's surface

In the triangle (PSG) in Figure (Rayp3d), the following expression relates the horizontal traveling of the source and geophone rays with the offset:  

$$\stackrel{\longrightarrow}{SG} = \stackrel{\longrightarrow}{SP} + \stackrel{\longrightarrow}{PG} .$$ (17)

The vectors $\stackrel{\longrightarrow}{PS}$ and $\stackrel{\longrightarrow}{PG}$ are related to the table that gives the lateral distance traveled by the ray as a function of the ray parameter and the travel time, $\xi(p,t)$. The projection on the earth's surface of a ray path (of parameter p and traveltime t) is the 2-D vector $\xi(p,t)\frac{{\bf p}}{p}$. Thus, relation (17) yields  

$${\bf x}_g - {\bf x}_s = \xi(p_g,t_g)\frac{{\bf p}_g}{p_g} - \xi(p_s,t_s)\frac{{\bf p}_s}{p_s} ,$$ (18)

where ${\bf x}_s$ and ${\bf x}_g$ are the 2-D vectors of the source and geophone coordinates.

A second relation expresses the coordinates of the point of emergence of the zero-offset ray, E. In the triangle (PEM), we have :  

$$\stackrel{\longrightarrow}{ME} = \stackrel{\longrightarrow}{MP} + \stackrel{\longrightarrow}{PE} .$$ (19)

Again, the vector $\stackrel{\longrightarrow}{PE}$ is related to the table of lateral distances, $\xi(p,t)$, yielding  

$${\bf x}_0 = - \frac{1}{2} \left( \xi(p_s,t_s)\frac{{\bf p}_... ...{{\bf p}_g}{p_g} \right) + \xi(p_0,t_0)\frac{{\bf p}_0}{p_0} ,$$ (20)

where ${\bf x}_0$ is the 2-D vector of the emergence point coordinates.

 

Rayp3d
Figure 22 Three-dimensional view of the source, receiver, and reflection points for 3-D v(z) DMO. The dashed lines represent the ray paths in the earth, and the bold solid lines represent the distances on the surface of the earth.

Time relationships in the earth's interior

The following obvious equation relates the traveltimes along the ray paths to source and geophone, t_s and t_g, to the given total traveltime t_sg:

 

t_s + t_g = t_sg . (21)

We can also state two more equations to account for the fact that the three ray paths, namely the zero-offset, source, and geophone ray paths, end at the same point in depth (R):  

$$\tau(p_0,t_0) = \tau(p_s,t_s) = \tau(p_g,t_g) ,$$ (22)

where $\tau(p,t)$ is the table of vertical traveltimes.

The relationship between the ray paths

At the reflector point R, the ray must obey Snell's law. In other words, the angle of incidence must equal the angle of reflection. In terms of ray paths, ${\bf r}_0$ fixes the position of the reflection plane, and ${\bf r}_s$ and ${\bf r}_g$ account for the angles of incidence and reflection. Since ${\bf r}_s$ and ${\bf r}_g$ have equal lengths 1/v, Snell's law can be expressed in the following vectorial form:  

$${\bf r}_0 = \lambda ({\bf r}_s + {\bf r}_g) ,$$ (23)

where $\lambda$ is a coefficient of proportionality. The first two coordinates are the ray parameter components, which have the following relation:  

$${\bf p}_0 = \lambda ({\bf p}_s + {\bf p}_g) .$$ (24)

The third coordinate of the ray path vectors can be expressed as a function of the inclination angle and the vector length, 1/v, as follows:  

$$r_z = \frac{\cos \theta(p,t)}{v} ,$$ (25)

where $\theta(p,t)$ is the table of angles along a ray of parameter p and time t. Then, after simplification of the 1/v factors (the velocity is the same for all rays at the reflector point), the third equation of relation (23) becomes  

$$\cos \theta(p_0,t_0) = \lambda ( \cos \theta(p_s,t_s) + \cos \theta(p_g,t_g) ) .$$ (26)

Substituting equation (26) into relation (24), we can eliminate the proportionality factor $\lambda$ to obtain the following relation:  

$${\bf p}_0 \left( \cos \theta(p_s,t_s) + \cos \theta(p_g,t_g ) \right) = ({\bf p}_s + {\bf p}_g) \cos \theta(p_0,t_0) .$$ (27)

Sets of equations and unknowns Our system of equations is constituted by collecting relations (18), (20), (21), (22), and (27):  

$$\left\{ \begin{array} {ccc} \xi(p_g,t_g)\frac{{\bf p}_g}{p... ... p}_s + {\bf p}_g) \cos \theta(p_0,t_0) \end{array} \right. .$$ (28)

Because the ${\bf x}$ and ${\bf p}$ vectors are two-dimensional, the first, second, and sixth relations of system (28) give $3\times2$ equations, yielding a set of nine equations.

The unknowns are ${\bf p}_s$, ${\bf p}_g$, t_s, t_g, t₀, and ${\bf x}_0$. As described in the preamble, the ${\bf p}$ vectors have only two unknown parameters, p_x and p_y (or p and $\phi$). Similarly, ${\bf x}_0$ is a two-dimensional unknown vector. Therefore, system (28) relates nine unknowns with the known parameters ${\bf p}_0$, t_sg, ${\bf x}_s$, ${\bf x}_g$, and the ray tracing tables $\xi(p,t)$, $\tau(p,t)$, and $\theta(p,t)$.

Solving the system

The second relation of system (28) can be isolated because it simply adds two relations and two unknowns (the coordinates of ${\bf x}_0$). Similarly, the third equation can be substituted into the others by replacing t_s (or t_g) with t_sg - t_g (or t_sg - t_s).

Thus, we obtain a reduced system of four relations (six equations):  

$$\left\{ \begin{array} {ccc} \xi(p_g, t_g)\frac{{\bf p}_g}{... ... p}_s + {\bf p}_g) \cos \theta(p_0,t_0) \end{array} \right. .$$ (29)

and two additional relations to compute the remaining unknowns:  

$${\bf x}_0 = \xi(p_0,t_0)\frac{{\bf p}_0}{p_0} - \frac{1}{2}... ...c{{\bf p}_s}{p_s} + \xi(p_g,t_g)\frac{{\bf p}_g}{p_g} \right)$$ (30)

 

t_s = t_sg - t_g . (31)

The six unknowns of system (29) are t₀, t_g (or t_s), ${\bf p}_s$ [including the two unknowns (p_xs,p_ys) or $(p_s,\phi_s)$], and ${\bf p}_g$ [(p_xg,p_yg) or $(p_g,\phi_g)$]. The parameters are t_sg, ${\bf p}_0$[including the two parameters (p_x0,p_y0) or $(p_0,\phi_0)$], ${\bf x}_s$, and ${\bf x}_g$. This system, like its 2-D counterpart, can be solved using the Newton-Raphson algorithm ().

Agreement with the equations in two dimensions

The 2-D system derived by Artley is a set of four equations with four unknowns. It can be retrieved by replacing the vectorial ray parameters with their scalar form (since the vectors are now in a common plane). The first relation of system (28) then becomes  

$$\xi(p_g,t_g) - \xi(p_s,t_{sg}-t_g) = 2h ,$$ (32)

where h is the half offset. The fourth relation becomes  

$$\frac{\sin \theta(p_0,t_0)}{v} \left( \cos \theta(p_s,t_s) ... ... + \frac{\sin \theta(p_g,t_g)}{v} \right) \cos \theta(p_0,t_0)$$ (33)

and can be simplified to  

$$\sin( \theta(p_0,t_0) - \theta(p_s,t_s) ) = \sin( \theta(p_g,t_g) - \theta(p_0,t_0) ) .$$ (34)

Finally, we obtain the bisection condition as shown by Artley , as follows:  

$$\theta(p_0,t_0) = \frac{1}{2}\left[ \theta(p_s,t_s) + \theta(p_g,t_g) \right] .$$ (35)

Conclusion We derived a set of equations that constrains the travel path of a ray in a three-dimensional earth model with a depth-variable velocity. This system can be solved for the zero-offset traveltime and the emergence point of the ray, yielding an efficient method for dip moveout processing. A successful implementation of the equations have been realized by Artley at the Colorado School of Mines, showing the expected saddle-like impulse response.

CONCLUSIONS The dip moveout correction requires important compromises between speed and precision in two- and three-dimensional models.

The conclusions to part I (page ) propose ways of increasing the speed. The constant-velocity assumption allows the dip moveout process to be fast in both two- and three-dimensional model geometries. Some additional attributes of the operator, such as weighting and anti-aliasing schemes, reduce the amplitude artifacts in a computationally efficient manner. Moreover, parallel computing can speed up the process even in the case of three-dimensional land data where the acquisition geometry is irregular. A parallel implementation of integral dip moveout in time slices proves efficient for irregular azimuthal distribution of the data.

The second part of this paper proposes two ways of improving precision. The first is to apply a proper weight along the operator, which is achieved at almost no extra computational cost. The second way is to consider depth-variable velocity. The dip moveout correction then becomes computationally expensive. However, an approximation valid for gently dipping reflectors allows the variable-velocity process to be almost as fast as the constant-velocity process. Aside from the considerations of speed, the method of three-dimensional dip moveout in depth-variable velocity developed in the last section may considerably improve the imaging of complex data structures.

Throughout my analysis of the dip moveout correction, I have tried to optimize both speed and precision. Unfortunately, the very fast methods prove inaccurate, whereas the very precise ones are slow. The approximate methods related to the squeezing of the dip moveout operator turn out to be fast and rather precise, even in three dimensions. This result suggests a method for three-dimensional dip moveout processing: the operator with limited cross-line extension that I introduced in the first part can be squeezed, yielding a first-order approximation of the theoretic ``saddle'' operator.

ACKNOWLEDGMENTS I would like to thank David Lumley for many useful discussions and for providing me with synthetic data and code. I am grateful to Mihai Popovici and Biondo Biondi for their efforts to transmit their knowledge of DMO and for their useful hints on parallel programming. I also thank Dave Nichols for his patience in teaching me the basis of parallel processing. Dimitri Bevc and David Lumley helped clarify for me the idea of anti-aliasing with triangle convolution. In addition, I would like to thank Craig Artley, Mihai Popovici, and Matthias Schwab for working with me on the derivation of the equations of 3-D v(z) DMO. Finally, I am grateful to Professor Jon Claerbout and Stew Levin for providing an excellent environment for my research, and to the Total oil company for funding it.

[SEP,SEGCON,SEGBKS,GEOTLE,EAEG,MISC,mybib]