## 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)
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)
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

## 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

### 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

### 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 $t0\u2212x0$ (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.