ABSTRACT

Time migration, as opposed to depth migration, suffers from two well-known shortcomings: (1) approximate equations are used for computing Green’s functions inside the imaging operator and (2) in case of lateral velocity variations, the transformation between the image-ray coordinates and the Cartesian coordinates is undefined in places where the image rays cross. We have found that the first limitation can be removed entirely by formulating time migration through wave propagation in image-ray coordinates. Our approach constructs a time-migrated image without relying on any kind of traveltime approximation by formulating an appropriate geometrically accurate acoustic wave equation in the time-migration domain. The advantage of this approach is that the propagation velocity in image-ray coordinates does not require expensive model building and can be approximated by quantities that are estimated in conventional time-domain processing. Synthetic and field data examples demonstrate the effectiveness of the proposed approach and show that the proposed imaging workflow leads to a significant uplift in terms of image quality and can bridge the gap between time and depth migrations. The image obtained by the proposed algorithm is correctly focused and mapped to depth coordinates; it is comparable to the image obtained by depth migration.

INTRODUCTION

Seismic imaging has the ultimate goal of creating an image of subsurface structures in depth (Gray et al., 2001; Bednar, 2005; Biondi, 2006; Etgen et al., 2009). However, it often falls short of this goal in practice because of the requirement of having an accurate velocity model (Glogovsky et al., 2008, 2009). That is possibly the main reason why time-domain imaging methods continue to play an important role in practical seismic data analysis.

Time migration takes a different route from depth migration by reducing the problem of velocity estimation to a parameter-picking problem. Each image point in a time-migrated image is associated with its own velocity parameter, which can be determined either by scanning different velocities (Yilmaz et al., 2001) or by wave extrapolation in the image-velocity space (Fomel, 2003). The computational advantage of time-domain imaging comes at the expense of two main flaws:

  1. 1)

    Time migration uses approximate Green’s functions (Zhang and Zhang, 1998) that typically rely on 1D velocity models and hyperbolic or slightly nonhyperbolic traveltime approximations.

  2. 2)

    In the case of lateral velocity variations, the transformation between the image-ray coordinates and Cartesian coordinates becomes distorted in places where the image rays cross. In such areas, there is no longer a one-to-one mapping between image-ray coordinates and Cartesian coordinates, and the coordinate transformation will also have a zero determinant (at the caustics of the image-ray field) (Hubral, 1977).

Image rays are seismic rays orthogonal to the surface of observation. These rays remain straight in the absence of lateral velocity variations but bend when they meet lateral heterogeneities. Cameron et al. (2007, 2008b, 2009) extend the image-ray theory to establish an exact theoretical connection between depth- and time-migration velocities and an inversion algorithm for converting the latter to the former. In the absence of lateral velocity variations, the time-to-depth conversion is accomplished by Dix inversion (Dix, 1955). As shown by Cameron et al. (2007) and Iversen and Tygel (2008), an additional correction is required when velocities vary laterally, that is, the correction related to the geometric spreading of image rays. Li and Fomel (2013, 2015) develop a robust algorithm for time-to-depth conversion including a geometric-spreading correction in the presence of lateral variations. Sripanich and Fomel (2018) develop a fast version of the time-to-depth conversion algorithm in the case of weak lateral variations.

In this paper, we propose to apply wave-equation imaging to create accurate seismic images in image-ray coordinates but without relying on Green’s function approximations and thus avoiding any inaccuracy issues associated with time migration (Fomel, 2013). We use the method of Sava and Fomel (2005) to define wave propagation in an alternative coordinate system and to connect it with the theory that relates time-migration velocities to velocity models in depth (Cameron et al., 2007). We show that, when the wave equation is transformed into the image-ray coordinate system, its coefficients are related to ideal time-migration velocities. Therefore, accurate one- or two-way wave-equation imaging can be accomplished by using information that is readily obtainable from conventional time-domain processing. This observation leads to a new imaging workflow, which provides a seamless transition from wave-equation time migration in the time domain to wave-equation time migration in the depth domain along with velocity model building in the depth domain for depth imaging. We test the proposed approach using simple synthetic and field data examples.

WAVE-EQUATION AND IMAGE RAYS

The phase information contained in seismic waves is governed by the eikonal equation, which takes the form (Chapman, 2004)
T·T=1V2(x),
(1)
where x is a point in space, T(x) is the traveltime from the source to x, and V(x) is the phase velocity. In an anisotropic medium, the phase velocity depends additionally on the direction of traveltime gradient T. For simplicity, we limit the following discussion to the isotropic case and 2D. Extensions to anisotropy and 3D are possible but would complicate the discussion. In 2D, x={x,z}, and the isotropic eikonal equation 1 can be written as
(Tx)2+(Tz)2=1V2(x,z).
(2)
The simplest wave equation that corresponds to eikonal equation 1 is the acoustic wave equation
2Wx2+2Wz2=1V2(x,z)2Wt2,
(3)
where W(x,z,t) is a wavefield propagating with velocity V. Equation 3 provides the basis for a variety of wave-equation migration algorithms, from one-way wave extrapolation in depth to two-way reverse time migration (RTM) (Biondi, 2006; Etgen et al., 2009).
Consider a family of image rays (Hubral, 1977), traced orthogonally to the surface. Image-ray coordinates x0 and t0 as functions of the Cartesian coordinates x and z satisfy the following system of partial differential equations (PDEs) (Cameron et al., 2007):
|x0|2=(x0x)2+(x0z)2=1J2(x,z),
(4)
x0·t0=x0xt0x+x0zt0z=0,
(5)
|t0|2=(t0x)2+(t0z)2=1V2(x,z),
(6)
with boundary conditions x0(x,0)=x and t0(x,0)=0. Equation 6 is the familiar eikonal equation 2. Equation 5 expresses the orthogonality of rays and wavefronts in an isotropic medium. Equation 4 defines the quantity J(x,z) as the geometric spreading of image rays.
Applying a change of variables from {x,z} to {x0,t0} transforms eikonal equation 2 to the coordinate system of image rays and leads to an elliptically anisotropic eikonal equation, which is
(Tx0)21J2+(Tt0)21V2=1V2.
(7)
The corresponding wave equation is given by
2Wx02Vd2+2Wt02=2Wt2,
(8)
where Vd=V/J. Equation 8 is a particular version of the more general Riemannian coordinate transformation analyzed by Sava and Fomel (2005). Additional versions of Riemannian coordinate transformations are analyzed by Shragge (2008) and Shragge and Shan (2008). Note that at the surface of observation z=0, the solution of equation 8 coincides geometrically with the solution of the original Cartesian equation 3.
The significance of equation 8 lies in the following fact established by Cameron et al. (2007). When time migration is performed using coordinates x0 and t0 and the conventional traveltime approximation based on the Taylor expansion of diffraction traveltime around the image ray, such as the classic hyperbolic approximation
t2(x)t02+(xx0)2V02(x0,t0),
(9)
the coefficient Vd appearing in equation 8 is related to time-migration velocity V0 appearing in equation 9. More specifically,
Vd2(x0,t0)=V2J2=[t0V02(x0,t0)]t0.
(10)
Following Cameron et al. (2008a), we call Vd the Dix velocity because it corresponds to the classic Dix inversion applied to the time-migration velocity (Dix, 1955). In the absence of lateral velocity variations (a V(z) medium), image rays do not spread; therefore, J=1, Vd=V and it corresponds to the true velocity in the medium, and V0 corresponds to the root-mean-square (rms) velocity. In the presence of lateral velocity variations, the conversion between time-migration coordinates and true Cartesian coordinates is more complicated (Bevc et al., 1995; Cameron et al., 2008a; Li and Fomel, 2015). However, to extrapolate seismic wavefields in image-ray coordinates using equation 8, it is sufficient to use an estimate of the Dix velocity Vd, which is readily available from conventional time-domain processing and equation 10.

WORKFLOW: WAVE-EQUATION TIME MIGRATION

We propose the following workflow for bridging the gap between time- and depth-domain imaging.

Step 1. Time-migration velocity analysis

Automatic velocity analysis in the process of time migration can be accomplished either by scanning a set of different velocities (Yilmaz et al., 2001) or by wave extrapolation in the image-velocity space (Fomel, 2003). It is also possible to estimate the time-migration velocity after the stack or from limited-offset data by separating and imaging the seismic diffractions (Harlan et al., 1984; Fomel et al., 2007; Burnett and Fomel, 2011; Dell and Gajewski, 2011; Coimbra et al., 2013; Decker et al., 2017). Alternatively, the time-migration velocity can be estimated directly from prestack data as a data attribute controlled by local slopes of reflection events (Fomel, 2007b; Cooke et al., 2009). Additional work on time-migration velocity analysis includes that by Gelius and Tygel (2015), Santos et al. (2015), and Glöckner et al. (2016).

Step 2. Dix conversion

Applying equation 10 for converting rms velocity V0(x0,t0) to Dix velocity Vd(x0,t0) in practice requires special care because of possible noisy measurements. Large variations in rms velocities produce rapid variations in interval velocities, which leads to unstable Dix inversion. It helps to formulate Dix conversion as a linear estimation problem and to use regularization for constraining its solution (Clapp et al., 1998; Valenciano et al., 2004; Fomel and Guitton, 2006; Fomel, 2007a). We use least squares and shaping regularization for stable conversion of rms velocities to interval velocities.

Step 3. Wave-equation time migration

Any of the available imaging techniques, such as Kirchhoff migration, one-way wave-equation migration, or two-way RTM, can be used to perform seismic imaging in image-ray coordinates using equation 8. Moreover, staying in this coordinate system allows migration velocity analysis and model building to be performed for estimating the Dix velocity Vd(x0,t0) instead of the usual seismic velocity V(x,z).

Step 4. Conversion from time to depth

Conversion from time to depth coordinates and from Vd(x0,t0) to V(x,z) is a nontrivial inverse problem. The problem involves not only a coordinate transformation (Hatton et al., 1981; Larner et al., 1981) but also a correction for the geometric spreading of image rays. As shown by Cameron et al. (2009), the problem can be reduced to solving an initial-value (Cauchy) problem for an elliptic PDE, which is a classic example of a mathematically ill-posed problem. To arrive at this formulation, let us transform the system of equations 46 into the image-ray coordinate system. The system of equations for inverse functions takes the form (Li and Fomel, 2015)
(xx0)2+(zx0)2=J2,
(11)
xx0xt0+zx0zt0=0,
(12)
(xt0)2+(zt0)2=V2,
(13)
with boundary conditions x(x0,0)=x0 and z(x0,0)=0.
From equation 12, it follows that
xt0=zx0zt0xx0.
(14)
Substituting this expression into equation 13 and using equation 11 leads to
xt0=Vd(x0,t0)zx0,
(15)
zt0=VJxx0=Vd(x0,t0)xx0,
(16)
where z/t0 and x/x0 are assumed to remain positive (image rays propagate down and do not cross). Finally, decoupling the system by using the equivalence of the second-order mixed derivatives produces the following system of two linear elliptic PDEs:
x0(Vdxx0)+t0(1Vdxt0)=0,
(17)
x0(Vdzx0)+t0(1Vdzt0)=0.
(18)
Additional initial conditions,
xt0|t0=0=0,
(19)
zt0|t0=0=Vd(x0,0),
(20)
specify that the image rays propagate down normal to the surface. If it were possible to solve systems 17 and 18 directly using only the initial conditions, the shape of the image rays could be determined and the true velocity could be estimated from equation 13. Unfortunately, this problem is mathematically ill-posed, which leads to numerical instability (Tikhonov and Arsenin, 1977). It can be approached, however, through regularization techniques (Cameron et al., 2009).
Li and Fomel (2015) develop robust time-to-depth conversion, which uses equations 46 in the Cartesian coordinate system and formulates time-to-depth conversion as a regularized least-squares optimization problem. Using linearization with respect to velocity perturbations, Sripanich and Fomel (2018) reformulate equations 17 and 18 for fast time-to-depth conversion appropriate for handling weak lateral variations. Weak lateral variation assumption is important because, in case of strong lateral variations, there is no longer a one-to-one mapping between image-ray coordinates and Cartesian coordinates, and the coordinate transformation will also have a zero determinant (at the caustics of the image-ray field). Using Sripanich and Fomel (2018), the squared Dix velocity converted to depth wd(x,z) is given as
wd(x,z)wdr(x,z)+(Δx0(x,z)×wdx0(x,z))+(Δt0(x,z)×wdt0(x,z)),
(21)
where wdr(x,z) denotes the wd(x0,t0) converted to depth based on the laterally homogeneous background assumption and the derivatives with respect to x0 and t0 are evaluated first in the original (x0,t0) coordinates followed by a similar conversion.

Step 5. Velocity model building

Finally, equipped with an estimate of the seismic velocity in Cartesian coordinates and a well-focused image, one could continue to refine the velocity model by using any of the conventional velocity estimation techniques (Robein, 2003; Jones, 2010).

EXAMPLES

Linear-gradient model

We start with a toy example of a linear-gradient model similar to the one used by Baina et al. (2002). The velocity in Figure 1a changes with a constant gradient tilted at 45°. In this model, the exact velocity is given by
v(x,z)=v0+gzz+gxx,
(22)
where v0=1 km/s, gz=0.15 1/s, and gx=0.15 1/s. Four reflectors with varying shapes are embedded in the model. Reflection data in Figure 1b are modeled using Kirchhoff modeling. The migration velocity squared wm and its Dix-inverted counterpart wd are given by the following expression (Li and Fomel, 2015):
wm(x0,t0)=((v0+gxx0)2t0(gcoth(gt0)gz))2,
(23)
wd(x0,t0)=((v0+gxx0)ggcosh(gt0)gzsinh(gt0))2.
(24)
The time-migration velocity computed using the analytical expression from Li and Fomel (2015) is shown in Figure 2a. The Kirchhoff time migration in Figure 2b fails to focus the image accurately because of strong lateral velocity variations. Figure 3a shows the Dix velocity that we further use to obtain the image by wave-equation time migration using RTM in image-ray coordinates (Dell et al., 2014) as shown in Figure 3b. The image is correctly focused but distorted because of image-ray bending. Bending image rays in Figure 4a correspond to the time-migration velocity shown in Figure 2a. The analytical solutions to time-to-depth conversion are (Li and Fomel, 2015)
x0(x,z)=x+(v0+gxx)2+gx2z2(v0+gxx)gx,
(25)
t0(x,z)=1garccosh[g2((v0+gxx)2+gx2z2+gzz)vgz2vgx2].
(26)
Figure 4b shows time-migration images converted to Cartesian coordinates. The image by wave-equation time migration is now well focused, correctly positioned in depth, and comparable in quality to the depth-migrated image in Figure 4c. Both images are created using low-rank RTM (Fomel et al., 2013).

Nankai field data example

To test the algorithm with field data, we first start with the Nankai data set (Forel et al., 2005). We preprocess the data set to correct for uneven bathymetry, ground-roll attenuation, and surface-consistent amplitude correction. To obtain the migration velocity, we use the Fowler (1984) dip moveout. The resultant migration velocity in Figure 5b is used to perform Kirchhoff time migration in Figure 5c of the stacked section in Figure 5a. Next, we convert the migration velocity to the Dix velocity in Figure 6a and subsequently use it for wave-equation time migration in Figure 6b. Magnified sections of the conventional time migration and wave-equation time migration (Figure 7) show that the image obtained by wave-equation time migration is correctly focused near the water column with the faults and subduction zone migrated to their true subsurface locations whereas time migration fails to focus the image accurately because of the strong lateral velocity variations. The image obtained by wave-equation time migration is still in the time coordinates. We transform it to the depth coordinates using the fast time-to-depth conversion algorithm (Sripanich and Fomel, 2018). Figure 8 shows Dix-inverted migration velocity squared wdr(x,z) and its gradients evaluated in the time-domain coordinates that are used to compute the interval velocity by the fast time-to-depth algorithm along with the image-ray coordinate system as shown in Figure 9a and 9b, respectively. Applying the time-to-depth conversion to the image obtained after the wave-equation time migration in Figure 10a and comparing it to the depth-migrated image in Figure 10b (obtained using the estimated interval velocity with the time-to depth conversion from the Dix velocity models), we see that the results are comparable and that the wave-equation time migration image is correctly focused and correctly positioned in depth.

Gulf of Mexico field data example

For the final example, we use a Gulf of Mexico field data set (Claerbout, 1995). In this data set, the maximum recording time is 4.0 s with the maximum offset of 3.48 km. The stacked section along with the picked migration velocity and Dix velocity is shown in Figure 11. We estimate the initial wdr(x,z) automatically using the method of velocity continuation (Fomel, 2003) followed by 1D Dix inversion to depth, which is similar to the workflow for time-to-depth conversion followed by Li and Fomel (2015) and Sripanich and Fomel (2018). Using the Dix velocity, we perform wave-equation time migration as shown in Figure 12b, which shows the improved delineation of faults compared to the conventional Kirchhoff time migration shown in Figure 12a. To convert the migrated section after wave-equation time migration to the depth domain, we use a fast time-to-depth conversion algorithm, similar to the previous example. The inputs to this model are shown in Figure 13. The output of the fast time-to-depth algorithm is the grid t0x0 (Figure 14a), which we use to map the wave-equation time migration results from time to depth coordinates (Figure 15a) and the estimated interval velocity in Figure 14b, which we use for depth migration in Figure 15b. We compare the final seismic image after a time-to-depth conversion process using wave-equation time migration with RTM. The results are comparable to Figure 15b, verifying the effectiveness of the proposed approach.

CONCLUSION

The proposed seismic imaging workflow is applicable to areas with mild lateral velocity variations. Mild variations are required to ensure that image rays do not cross, so that the mapping between image-ray coordinates and Cartesian coordinates remains one to one. We have outlined a theoretical foundation for the central step of this workflow (wave-equation time migration) and demonstrated its application using synthetic and field data examples. Wave-equation time migration produces an image in time-migration (image-ray) coordinates but without using moveout approximations or any other limitations associated with traditional prestack-time-migration algorithms. The proposed algorithm can be cost effective because it circumvents the need for multiple migrations that are required to update the velocity model in depth for the conventional algorithms during the velocity model building process. In complex velocity models, when the image-ray coordinate system breaks down because of crossing rays, the proposed workflow may not be directly applicable. However, it can be combined with redatuming (downward extrapolation) to allow the velocity model to be recursively estimated below the depth of the breakdown point.

ACKNOWLEDGMENTS

We thank the sponsors of Texas Consortium of Computational Seismology for financial support.

DATA AND MATERIALS AVAILABILITY

Data associated with this research are available and can be obtained by contacting the corresponding author.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.