Diffraction imaging and time-migration velocity analysis using oriented velocity continuation

We perform seismic diffraction imaging and time-migration velocity analysis by separating diffractions from specular reflections and decomposing them into slope components. We image the slope components using migration velocity extrapolation in time-space-slope coordinates. The extrapolation is described by a convection-type partial differential equation and implemented in a highly parallel manner in the Fourier domain. Synthetic and field data experiments show that the proposed algorithms are able to detect accurate timemigration velocities by measuring the flatness of diffraction events in slope gathers for singleand multiple-offset data.


INTRODUCTION
Seismic diffraction occurs when a seismic wave encounters a heterogeneity without a clearly defined tangent plane, such as an edge or tip, and the reflection part of the ray theory breaks down (Klem-Musatov, 1994). These divergent diffraction rays have similar behavior to a subsurface secondary source located at the heterogeneity (Keller, 1962). Analyzing diffraction moveout behavior in different domains can provide subsurface velocity information analogous to analyzing reflection moveout behavior from a surface source.
The fact that diffractions migrated with correct velocity collapse to points motivated Harlan et al. (1984) to propose the idea of separating diffractions from specular reflections and using diffraction focusing as a tool for velocity analysis. Separation of diffraction events from seismic data is a necessary step for velocity analysis because diffraction signals are typically significantly weaker than those of reflections (Klem-Musatov, 1994). Fomel et al. (2007) developed a constructive procedure for diffraction separation based on plane-wave destruction and diffraction focusing analysis based on velocity continuation and local kurtosis. The procedure was extended to 3-D azimuthally-anisotropic velocity analysis by Burnett and Fomel (2011). However, local kurtosis may not be an optimal measure for diffraction focusing because it requires smoothing or windowing in space, which reduces spatial velocity resolution through the smoothing window parameters.

TCCS-12 Oriented velocity continuation
A particularly convenient domain for separating diffractions and reflections and for analyzing migration velocities is dip-angle gathers Biondi and Symes, 2004;Landa et al., 2008;Reshef and Landa, 2009;Klokov and Fomel, 2012). In the dip-angle domain, specular reflections appear as hyperbolic events centered at the reflector dip and bending upwards, even when over or undermigrated, and diffractions appear flat when imaged at the location of the diffractor with the correct velocity (Reshef, 2007). Measuring flatness of diffraction events in dip-angle gathers, as opposed to flatness of reflection and diffraction events in reflection-angle gathers, provides an alternative constraint on seismic velocity (Reshef and Landa, 2009). Traditionally, dip-angle gathers are constructed with Kirchhoff migration (Fomel and Prucha, 1999;Xu et al., 2001;Cheng et al., 2011;Koren and Ravve, 2011;Bashkardin et al., 2012;Klokov and Fomel, 2013).
In this paper, we adopt an analogous method to the dip-angle approach used by Reshef and Landa (2009) to devise a constructive and highly parallel procedure for estimating velocities in time-domain processing using data decomposition in slope (Ghosh and Fomel, 2012) and velocity continuation in the midpoint-time-slope domain. By analogy with the "oriented wave equation" (Fomel, 2003a), we call this approach oriented velocity continuation (OVC ) and develop a fast spectral method for its implementation on common-offset data. This differs from the methods devised by Reshef and Landa (2009), which utilize a separate Kirchhoff-based angle prestack time or depth migration and calculation of travel time tables for each tested migration velocity. OVC uses a continuation approach where a single migration is used to determine an initial image in the midpoint-time-slope domain to which a velocity dependent phase shift is applied over the range of plausible migration velocities, enabling OVC to test a greater number of velocities at a lower computational cost.
Using a field-data experiment, we demonstrate the effectiveness of oriented velocity continuation in zero-offset diffraction imaging and velocity analysis, and using a synthetic model we observe that higher velocity resolution can be achieved when multiple offsets are included in the process.

ORIENTED VELOCITY CONTINUATION
Velocity continuation (Fomel, 2003c) is the imaginary process of a continuous transformation of seismic time-migrated images as they are propagated through different migration velocities. In the most general terms, the kinematics of velocity continuation can be described by an equation of the Hamilton-Jacobi type where τ (x, v) is the location of a time-migrated reflector with time-domain coordinates x = (x 1 , x 2 , t) imaged with spatially constant time-migration velocity v. The particular form of function F in equation (1) depends on the acquisition geometry of TCCS-12 Oriented velocity continuation the input data. For the case of common-offset 2D velocity continuation for data with half-offset h, and equation (1) corresponds to the characteristic equation of the image propagation process which describes a propagation of the time-migrated image I(t, x, v) in velocity v (Fomel, 2003c). Time-domain imaging can be performed effectively by extrapolating images in velocity and estimating velocity v m (t, x) of the best image (Larner and Beasley, 1987;Fomel, 2003b;Fomel and Landa, 2014).
As shown by Fomel (2003a), it is possible to extend the formulation of a wave propagation process from the usual time-and-space coordinates to the phase space consisting of time, space, and slope. Applying a similar approach to equation (1), we first employ the Hamilton-Jacobi theory (Courant and Hilbert, 1989;Evans, 2010) to write the corresponding system of ordinary differential equations for the characteristics (velocity rays), as follows: where p stands for ∇τ , the gradient or slope of time-migrated wavefield energy.
If the image I(t, x, v) is decomposed in slope components I(t, x, p, v) so that we can then look for an equation that would adequately describe a continuous transformation of I. To preserve the geometry of the transformation, it is sufficient to require that I transports along the characteristics described by equations (3-5). Applying partial derivatives and the chain rule, we arrive at the equation analogous to the Liouville equation (Engquist and Runborg, 2003): Equation (7) describes, in the most general form, the process of oriented velocity continuation, image propagation in velocity in the coordinates of time-space-slope. It is a linear first-order partial differential equation of convection type which operates in the phase space.

Common-offset oriented velocity continuation
To adopt the general theory described above to the case of common-offset 2D oriented velocity continuation, we can substitute equation (2) into (7), arriving at the equation which describes image propagation in the time-space-slope coordinates rather than the usual time-space coordinates. After this kind of extrapolation, regular images can be reconstructed by stacking over offset and slope.
Slope gathers, analogous to dip-angle gathers, can be extracted before stacking over slope by analyzing {t, p} panels for different image locations x and velocities v. Measuring flatness of diffraction events in these gathers provides a means for estimating migration velocity (Landa et al., 2008;Reshef and Landa, 2009).
For practical implementation, the formulation of oriented velocity continuation can be simplified by employing a stretch from the regular time coordinate to squared time σ = t 2 (Fomel, 2003b). According to this transformation, the Hamilton-Jacobi equation (1) which leads to the simpler form of the oriented equation where q corresponds to ∂σ ∂x , and the image is constructed in {σ, x, q, h} coordinates instead of {t, x, p, h} coordinates. Applying the Fourier transform, we can further transform equation (10) to Equation (11) has the analytical solution: where v 0 is a constant non-zero initial migration velocity.
Stacking over offset provides a slope-decomposed formulation for oriented velocity continuation: .

(13) Oriented velocity continuation
This derivation suggests the following algorithm for time-domain imaging using common-offset 2D oriented velocity continuation: 1. Start with initial time migration with a constant velocity v 0 to generate I(t, x, v 0 , h).
2. Apply vertical time stretch to transform from t to σ.
Note that this operation is parallel in ω and h.

Apply
Fourier transform from x to k to generateĨ(ω, k, q, v 0 , h). Note that this operation is parallel in q and h.
6. Apply the phase-shift filter from equation (12) to generateĨ(ω, k, q, v, h) for multiple values of v. Note that this operation is data-intensive but parallel in q, k, and h.
9. Apply inverse time stretch from σ to t.
10. Stack over q and extract the slice at time-migration velocity v m (t, x) to generate the final time-migrated image I(t, x, v m (t, x)).
In order to estimate the velocity v m (t, x), we apply the workflow described above to diffraction imaging and modify it as follows: • Before Step 1, we separate reflections and diffractions in the common-offset data using local plane-wave destruction (Fomel, 2002;Fomel et al., 2007;Decker et al., 2013).
• After Step 9, we analyze slope gathersÎ(t, x, q, v) and automatically pick the velocity v m (t, x) that corresponds to the maximum flatness (semblance) over q using the picking algorithm described by Fomel (2009b). This approach follows the principle of flatness of diffraction events in slope gathers (Landa et al., 2008;Reshef and Landa, 2009;Klokov and Fomel, 2012).
The computational cost associated with determining velocity using oriented velocity continuation is linear with the number of time samples, spatial samples, offsets, velocities, and slopes considered. It is parallel in spatial samples, offsets, velocity, and slope. The cost may then be considered as O NtNxN h NvNp Nc where N c is the number of cores available.

SLOPE DECOMPOSITION
In order to perform the initial slope decomposition (Step 4 in the algorithm above), we adopt the method of Ghosh and Fomel (2012). The idea of slope decomposition was discussed previously by Ottolini (1983) and implemented using the local-slant stack transform (Ventosa et al., 2012). The slope-decomposition algorithm suggested by Ghosh and Fomel (2012) is based on the time-frequency decomposition of Liu and Fomel (2013). Namely, at each frequency ω, we apply regularized non-stationary regression (Fomel, 2009a) to transform from space x to space-slope x-q domain. The non-stationary regression amounts to finding complex coefficients A n (ω, x) in the decomposition where D(ω, x) is the image slice, and D n (ω, x) is its slope component corresponding to slope q n : Equation (14) is the discrete analog of equation (6). Similarly to the time-frequency decomposition proposed by Liu and Fomel (2013), shaping regularization is used to control the variability of A n coefficients and to accelerate the algorithm.

EXAMPLES Toy Model Example
We first illustrate the concept of oriented velocity continuation using a simple toy model (Landa et al., 2008;Klokov and Fomel, 2013) with constant 1.0 km s velocity containing one dipping and one flat reflector and a single diffractor centered at 0.5 km (Figure 1a). Figure 1b shows data warped to squared time.
Warped data are decomposed into their constituent slope components and initially under-migrated using v 0 = 0.5 km s . The initial slope decomposed image is shown in Figure 2, which illustrates a slope gather centered above the diffractor on the right panel, and a partial image containing energy with the slope of the top dipping reflector on the front panel. The partial image contains energy of the top dipping reflector, which has the selected slope, diffraction energy with that slope, and a small portion of the energy from the flat bottom reflector. Stacking over all constituent slopes provides an image (top left panel of Figure 4). TCCS-12 Oriented velocity continuation Figure 2: Slope decomposed initial migration for toy model data using v 0 = 0.5 km s . The side pane shows the slope decomposed data centered at 0.5 km, the diffractor, and the front pane shows a partial image for data containing energy with a slope of 2.998 km s 2 , that of the top reflector.
Examining the slope gather of the initial migration in Figure 3a, the three panels contain points of energy corresponding, from top to bottom, to the top reflector, the diffractor, and the bottom reflector. The energy of each reflector is contained at the same lateral position in the three panels due to the constant slope of the reflector, although the vertical position of the top dipping reflector changes through the panels because the reflector dips downward to the right.
The diffraction has a hyperbolic moveout rather than a constant slope, so energy appears at different slopes in different slope gathers, with zero slope in the gather centered over the diffractor at 0.5 km. As we propagate data through velocity, this pattern holds: the reflection energy is stationary at its slope location for all gathers and diffraction energy has zero slope when viewed in the gather above the diffractor and non-zero slope for other gathers.
The initially migrated image is propagated to the higher time migration velocity of 0.75 km s using oriented velocity continuation, and slope gathers are illustrated in Figure 3b. Reflection energy now bends upward about the stationary point of each reflection in "smiles" that become more accentuated with larger migration velocities. The diffraction event bends upward as well, but this only holds for the current case of under-migration. Figure 3c shows the image propagated to the correct migration velocity. Diffrac-TCCS-12 Oriented velocity continuation tion energy is planar in all three gathers and flat in the middle gather centered above the diffractor. Stacking the energy in each gather over slope collapses the flat diffraction energy in the central gather to a point at the location of the diffractor. Sloping diffraction energy in the right and left panels cancels out when summed over slope. This flatness is essential to using oriented velocity continuation as a tool for determining the correct migration velocity.
When data are further propagated to over-migration with velocity of 1.25 km s in Figure 3d, the diffraction event bows downward in a "frown" juxtaposed against the upward bending reflection "smiles".
Stacking over slope provides images for these four velocities (Figure 4). In these images, the diffraction event incrementally evolves from having a hyperbolic downward character in the top under-migrated row, to collapsing to a point in the bottom left panel with the correct migration velocity, to bowing hyperbolically upward in the bottom right over-migrated image.
The changing geometry of diffraction energy in the gathers can be harnessed to determine the proper migration velocity (Landa et al., 2008;Reshef and Landa, 2009). The correct migration velocity will be the one that maximizes the "flatness" of slope decomposed diffraction events, as measured by coherence or another appropriate metric. Selecting 1.0 km s , the velocity that produces flat diffraction energy, provides us with a properly migrated image ( Figure 5).
To properly estimate migration velocity using diffraction flatness, reflection events must first be filtered out from the diffraction data, or else the contribution of reflection "smiles" may dominate and bias the flatness measure.

Synthetic Example
To test oriented velocity continuation and its velocity resolution we generate a synthetic dataset using a model with a constant velocity gradient beginning with a 2.0 km s surface velocity. Diffractors are created as reflectivity spikes within the model with random spatial and magnitude distributions. Kirchhoff forward modeling is used to generate 24 offsets with a 50 m interval.
Zero-offset data are shown in Figure 6a, and 1.0 km common-offset data appear in Figure 6b. A time shift between the data is noticeable.
Both zero and common-offset data are warped to squared time, slope decomposed, and migrated with a 2.0 km s initial velocity. The initially migrated slope decomposed images are propagated through a range of plausible migration velocities using oriented velocity continuation. Common-offset partial images are then stacked over offset for each continuation velocity. the correct velocity for the diffractor at 1.4 s located directly underneath the midpoint of this gather. When a different velocity is used, the shape of the event deviates from planar. We perform velocity analysis by testing the semblance, or flatness, of diffraction events in slope gathers over the range of velocities. Because velocity does not vary laterally in our synthetic model, we average the semblance across midpoints to generate semblance panels for the zero-offset and common-offset cases (Figures 7a  and 7b respectively).
As seen from the slope gathers (Figure 9), for the zero-offset case, there is a stationary point corresponding to diffraction energy with zero slope which does not shift vertically under velocity perturbations. For the 1.15 km common-offset case, perturbing velocity changes the slope decomposed diffraction shape and shifts it vertically. Slope gathers with incorrect velocities (Figure 9b) are time shifted with respect to those generated for the zero-offset case (Figure 9a). When the correct migration velocity is used, horizontal common-offset diffraction energy appears at the same time as for the zero-offset case. The vertical shift of incorrectly migrated common-offset data leads to a sharper change in estimated flatness values while converging on the correct migration velocity and therefore improves velocity resolution. Therefore, commonoffset semblance panel appears to have higher spatial and vertical resolution than the zero-offset case. This higher spatial resolution can be attributed to the improved illumination of scattering objects with the full range of offsets.
Final images for zero-offset and stacked common-offset cases using migration velocities estimated from the semblance panels are shown in Figure 8. Differences between the two images are too small to easily detect in this example. However, due to the higher velocity resolution visible in the semblance panel resulting from the consideration of multiple offsets, Figure 7b, we expect the stacked common-offset image to be better resolved and less prone to noise than zero-offset case when applied to field datasets.

Field Data Example
We demonstrate an application of zero-offset (post-stack) oriented velocity continuation on a deep water 2D line acquired to image the Nankai Trough subduction zone. Data acquisition parameters as well as processing results can be found in Moore et al. (1990), where the line is referred to as NT62-8. Structural interpretation can be found in Moore and Shipley (1993). Here, we consider a fragment of the line (CMPs 900-1301) used previously by Forel et al. (2005).
Conventional velocity analysis resolution suffers in this dataset from the limitations imposed by the depth of a seabed in the area (average of ≈ 4.5 km) and a relatively short 2 km streamer length. For deep water datasets diffractions may exhibit better illumination than reflections because diffraction aperture is not restricted to the recording array length, enabling them to provide a potentially more detailed velocity distribution. This behavior makes OVC migration velocity analysis appealing.
The DMO stacked section considered in this study is shown in Figure 10. Diffractions are extracted via plane-wave destruction (Figure 11), warped to squared time, and decomposed into slope. Figure 12 shows slope decomposed data warped back to regular time for ease of comparison with slope decomposed images appearing later. Next, we take the decomposed data through oriented velocity continuation over a range of sixty constant migration velocities beginning with v 0 = 1.4 km s using a 20 m s step. Diffraction events bend upward in the slope gather centered above x = 4100.06 m with the minimum tested migration velocity (Figure 14a), indicating under-migration. Diffraction events in the slope gather centered above the same location with the maximum tested migration velocity (Figure 14b) bend downward, indicating over-migration.
Gather semblance is calculated for each continuation velocity, and migration velocity is automatically picked by attempting to maximize semblance for plausible velocity values at each CMP location. Semblance panels with superimposed picks are shown in Figure 15. Anomalies corresponding to higher velocities than the picked trend may correspond to reflections with high curvature, like the one located between x = 3000 and 4000 m between t = 6.0 and 6.5 s in Figures 10 and 11. Highly curved reflections have a similar behavior to diffractions in response to migration velocity perturbation (Sava et al., 2005), but focus at a higher velocity than the correct one. Intersection of over-migrated reflection and diffraction tails from the rugose seabed, some of which are out of plane, leads to diffraction-like events, another cause of false semblance highs. These are visible in the three semblance panels of Figure 15 above 6.4 s. Low velocity semblance anomalies corresponding to the flattening of out of plane diffractions are also visible in the middle interval of the semblance panels, particularly near t ≈ 6.8 s in the central and right panels centered above x = 5000 and 6500 m.
Combining the semblance velocity picks from each CMP provides a time-migration velocity field, shown in Figure 16. As noted above, several anomalously low velocity TCCS-12 Oriented velocity continuation zones exist in the picked field, primarily between t = 6.5 and 7 s where the attempted flattening of out of plane diffractions leads to a low picked velocity.
Gathers corresponding to the picked velocity are selected. Examining a slope gather from x = 4100.06 m generated using the picked migration velocity (Figure 14c), diffraction events now appear flat, particularly the one located near t = 6.4 s , indicating that they have been correctly migrated.
Stacking gathers generated from the the picked velocity over slope provides the diffraction image in Figure 17. We apply oriented velocity continuation to the DMO stacked data from Figure 10 and stack over gathers selected with the appropriate velocity to generate the image of reflections and diffractions in Figure 18. Both images highlight fault surfaces. Finer discontinuities, such as those associated with the rough surface of the subducting plate crust, located near t ≈ 7.5 s (Moore and Shipley, 1993), are more prominent on the diffraction image and tend to be well focused, supporting the accuracy of the picked velocity.

CONCLUSIONS
We have developed and demonstrated a highly parallel and constructive procedure for time-domain velocity estimation. The method operates by decomposing data by Oriented velocity continuation  slope and propagating slope components in velocity in the midpoint-time-slope domain. Semblance in slope gathers is used as a measure for selecting velocities that correspond to correctly migrated flat diffraction events, which applies even for singleoffset data. This semblance measure is observed to achieve higher resolution when multiple offsets are considered in the oriented velocity continuation process. If data with multiple offsets are available, oriented velocity continuation could be used to better constrain the velocity model when used in conjunction with traditional reflection moveout analysis. Chosen velocities can be used to generate both diffraction and reflection images. The powerful ability of oriented velocity continuation to operate with only zero-offset data enables accurate migration velocity analysis in situations where only limited-offset data are available.
Flatness of slope-decomposed diffraction events is more responsive to velocity perturbation than diffraction focusing, because it does not require smoothing or windowing in space. Therefore, the proposed method has the potential for diffraction velocity estimation with superior resolution when compared to methods based on diffraction focusing.
Oriented Velocity Continuation is formulated as a type of time-migration, and is thus subject to constraints relating to image distortion from horizontal velocity changes in the subsurface. The presence of strong lateral velocity variations may alter diffraction moveout, making event slope change with azimuth. In such a case,  OVC would be unable to completely flatten the diffraction signal in slope gathers and locally determine the correct migration velocity.
This method can be extended to three dimensions using data decomposition by azimuth and inclination for each image point.
Operating on three dimensional data should improve velocity resolution by overcoming out of plane artifacts in the seismic image. Although extension to 3D adds the expense of additional spatial and slope dimensions, the Fourier-domain computation would also be parallel in these new dimensions, making the operation feasible in practice using computer clusters.