Wave-equation time migrationWave-equation time migration | Geophysics

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

The phase information contained in seismic waves is governed by the eikonal equation, which takes the form (Chapman, 2004)

\nabla T \cdot \nabla T = \frac{1}{V^{2} (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

\nabla 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 can be written as

(\frac{\partial T}{\partial x})^{2} + (\frac{\partial T}{\partial z})^{2} = \frac{1}{V^{2} (x, z)} .

(2)

The simplest wave equation that corresponds to eikonal equation is the acoustic wave equation

\frac{\partial^{2} W}{\partial x^{2}} + \frac{\partial^{2} W}{\partial z^{2}} = \frac{1}{V^{2} (x, z)} \frac{\partial^{2} W}{\partial t^{2}},

(3)

where

W (x, z, t)

is a wavefield propagating with velocity

V

⁠. Equation 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

x_{0}

and

t_{0}

as functions of the Cartesian coordinates

x

and

z

satisfy the following system of partial differential equations (PDEs) (Cameron et al., 2007):

| \nabla x_{0} |^{2} = (\frac{\partial x_{0}}{\partial x})^{2} + (\frac{\partial x_{0}}{\partial z})^{2} = \frac{1}{J^{2} (x, z)},

(4)

\nabla x_{0} \cdot \nabla t_{0} = \frac{\partial x_{0}}{\partial x} \frac{\partial t_{0}}{\partial x} + \frac{\partial x_{0}}{\partial z} \frac{\partial t_{0}}{\partial z} = 0,

(5)

| \nabla t_{0} |^{2} = (\frac{\partial t_{0}}{\partial x})^{2} + (\frac{\partial t_{0}}{\partial z})^{2} = \frac{1}{V^{2} (x, z)},

(6)

with boundary conditions

x_{0} (x, 0) = x

and

t_{0} (x, 0) = 0

⁠. Equation is the familiar eikonal equation . Equation expresses the orthogonality of rays and wavefronts in an isotropic medium. Equation defines the quantity

J (x, z)

as the geometric spreading of image rays.

Applying a change of variables from

{x, z}

{x_{0}, t_{0}}

transforms eikonal equation to the coordinate system of image rays and leads to an elliptically anisotropic eikonal equation, which is

(\frac{\partial T}{\partial x_{0}})^{2} \frac{1}{J^{2}} + (\frac{\partial T}{\partial t_{0}})^{2} \frac{1}{V^{2}} = \frac{1}{V^{2}} .

(7)

The corresponding wave equation is given by

\frac{\partial^{2} W}{\partial x_{0}^{2}} V_{d}^{2} + \frac{\partial^{2} W}{\partial t_{0}^{2}} = \frac{\partial^{2} W}{\partial t^{2}},

(8)

where

V_{d} = V / J

⁠. Equation 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 coincides geometrically with the solution of the original Cartesian equation .

The significance of equation lies in the following fact established by Cameron et al. (2007). When time migration is performed using coordinates

x_{0}

and

t_{0}

and the conventional traveltime approximation based on the Taylor expansion of diffraction traveltime around the image ray, such as the classic hyperbolic approximation

t^{2} (x) \approx t_{0}^{2} + \frac{(x - x_{0})^{2}}{V_{0}^{2} (x_{0}, t_{0})},

(9)

the coefficient

V_{d}

appearing in equation is related to time-migration velocity

V_{0}

appearing in equation . More specifically,

V_{d}^{2} (x_{0}, t_{0}) = \frac{V^{2}}{J^{2}} = \frac{\partial [t_{0} V_{0}^{2} (x_{0}, t_{0})]}{\partial t_{0}} .

(10)

Following Cameron et al. (2008a), we call

V_{d}

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

⁠,

V_{d} = V

and it corresponds to the true velocity in the medium, and

V_{0}

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 , it is sufficient to use an estimate of the Dix velocity

V_{d}

⁠, which is readily available from conventional time-domain processing and equation .

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 $V_{0} (x_{0}, t_{0})$ to Dix velocity $V_{d} (x_{0}, t_{0})$ 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 $V_{d} (x_{0}, t_{0})$ instead of the usual seismic velocity $V (x, z)$ ⁠.

Step 4. Conversion from time to depth

Conversion from time to depth coordinates and from

V_{d} (x_{0}, t_{0})

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 – into the image-ray coordinate system. The system of equations for inverse functions takes the form (Li and Fomel, 2015)

(\frac{\partial x}{\partial x_{0}})^{2} + (\frac{\partial z}{\partial x_{0}})^{2} = J^{2},

(11)

\frac{\partial x}{\partial x_{0}} \frac{\partial x}{\partial t_{0}} + \frac{\partial z}{\partial x_{0}} \frac{\partial z}{\partial t_{0}} = 0,

(12)

(\frac{\partial x}{\partial t_{0}})^{2} + (\frac{\partial z}{\partial t_{0}})^{2} = V^{2},

(13)

with boundary conditions

x (x_{0}, 0) = x_{0}

and

z (x_{0}, 0) = 0

⁠.

From equation , it follows that

\frac{\partial x}{\partial t_{0}} = - \frac{\frac{\partial z}{\partial x_{0}} \frac{\partial z}{\partial t_{0}}}{\frac{\partial x}{\partial x_{0}}} .

(14)

Substituting this expression into equation and using equation leads to

\frac{\partial x}{\partial t_{0}} = V_{d} (x_{0}, t_{0}) \frac{\partial z}{\partial x_{0}},

(15)

\frac{\partial z}{\partial t_{0}} = \frac{V}{J} \frac{\partial x}{\partial x_{0}} = V_{d} (x_{0}, t_{0}) \frac{\partial x}{\partial x_{0}},

(16)

where

\partial z / \partial t_{0}

and

\partial x / \partial x_{0}

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:

\frac{\partial}{\partial x_{0}} (V_{d} \frac{\partial x}{\partial x_{0}}) + \frac{\partial}{\partial t_{0}} (\frac{1}{V_{d}} \frac{\partial x}{\partial t_{0}}) = 0,

(17)

\frac{\partial}{\partial x_{0}} (V_{d} \frac{\partial z}{\partial x_{0}}) + \frac{\partial}{\partial t_{0}} (\frac{1}{V_{d}} \frac{\partial z}{\partial t_{0}}) = 0 .

(18)

Additional initial conditions,

\frac{\partial x}{\partial t_{0}} |_{t_{0} = 0} = 0,

(19)

\frac{\partial z}{\partial t_{0}} |_{t_{0} = 0} = V_{d} (x_{0}, 0),

(20)

specify that the image rays propagate down normal to the surface. If it were possible to solve systems and directly using only the initial conditions, the shape of the image rays could be determined and the true velocity could be estimated from equation . 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 – 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 and 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

w_{d} (x, z)

is given as

w_{d} (x, z) \approx w_{d r} (x, z) + (Δ x_{0} (x, z) \times \frac{\partial w_{d}}{\partial x_{0}} (x, z)) + (Δ t_{0} (x, z) \times \frac{\partial w_{d}}{\partial t_{0}} (x, z)),

(21)

where

w_{d r} (x, z)

denotes the

w_{d} (x_{0}, t_{0})

converted to depth based on the laterally homogeneous background assumption and the derivatives with respect to

x_{0}

and

t_{0}

are evaluated first in the original

(x_{0}, t_{0})

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 ¹ changes with a constant gradient tilted at 45°. In this model, the exact velocity is given by

v (x, z) = v_{0} + g_{z} z + g_{x} x,

(22)

where

v_{0} = 1

km/s,

g_{z} = 0.15

1/s, and

g_{x} = 0.15

1/s. Four reflectors with varying shapes are embedded in the model. Reflection data in Figure ¹ are modeled using Kirchhoff modeling. The migration velocity squared

w_{m}

and its Dix-inverted counterpart

w_{d}

are given by the following expression (Li and Fomel, 2015):

w_{m} (x_{0}, t_{0}) = (\frac{(v_{0} + g_{x} x_{0})^{2}}{t_{0} (g \coth (g t_{0}) - g_{z})})^{2},

(23)

w_{d} (x_{0}, t_{0}) = (\frac{(v_{0} + g_{x} x_{0}) g}{g \cosh (g t_{0}) - g_{z} \sinh (g t_{0})})^{2} .

(24)

The time-migration velocity computed using the analytical expression from Li and Fomel (2015) is shown in Figure ². The Kirchhoff time migration in Figure ² fails to focus the image accurately because of strong lateral velocity variations. Figure ³ 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 ³. The image is correctly focused but distorted because of image-ray bending. Bending image rays in Figure ⁴ correspond to the time-migration velocity shown in Figure ². The analytical solutions to time-to-depth conversion are (Li and Fomel, 2015)

x_{0} (x, z) = x + \frac{\sqrt{(v_{0} + g_{x} x)^{2} + g_{x}^{2} z^{2}} - (v_{0} + g_{x} x)}{g_{x}},

(25)

t_{0} (x, z) = \frac{1}{g} arccosh [\frac{g^{2} (\sqrt{(v_{0} + g_{x} x)^{2} + g_{x}^{2} z^{2}} + g_{z} z) - v g_{z}^{2}}{v g_{x}^{2}}] .

(26)

Figure ⁴ 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 ⁴. 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 ⁵ is used to perform Kirchhoff time migration in Figure ⁵ of the stacked section in Figure ⁵. Next, we convert the migration velocity to the Dix velocity in Figure ⁶ and subsequently use it for wave-equation time migration in Figure ⁶. Magnified sections of the conventional time migration and wave-equation time migration (Figure ⁷) 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 ⁸ shows Dix-inverted migration velocity squared $w_{d r} (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 ⁹ and ⁹, respectively. Applying the time-to-depth conversion to the image obtained after the wave-equation time migration in Figure ¹⁰ and comparing it to the depth-migrated image in Figure ¹⁰ (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 ¹¹. We estimate the initial $w_{d r} (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 ¹², which shows the improved delineation of faults compared to the conventional Kirchhoff time migration shown in Figure ¹². 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 ¹³. The output of the fast time-to-depth algorithm is the grid $t_{0} - x_{0}$ (Figure ¹⁴), which we use to map the wave-equation time migration results from time to depth coordinates (Figure ¹⁵) and the estimated interval velocity in Figure ¹⁴, which we use for depth migration in Figure ¹⁵. 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 ¹⁵, 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.