# Azimuthally Anisotropic 3D Velocity Continuation

Authors: William Burnett; Sergey Fomel
Journal: International Journal of Geophysics
 (2011)

---

## Abstract

We extend time-domain velocity continuation to the zero-offset 3D azimuthally anisotropic
case. Velocity continuation describes how a seismic image changes given a change in migration
velocity. This description turns out to be of a wave propagation process, in which images
change along a velocity axis. In the anisotropic case, the velocity model is multiparameter.
Therefore, anisotropic image propagation is multidimensional. We use a three-parameter
slowness model, which is related to azimuthal variations in velocity, as well as their principal
directions. This information is useful for fracture and reservoir characterization from
seismic data. We provide synthetic diffraction imaging examples to illustrate the concept
and potential applications of azimuthal velocity continuation and to analyze the impulse
response of the 3D velocity continuation operator.


---

## Body

## 1. Introduction

Velocity continuation [1, 2] provides a framework for describing how a seismic image changes given a change in the migration velocity model. Similar in concept to residual migration [3], and cascaded migrations [4], velocity continuation is a continuous formulation of the same process. Velocity continuation has found applications in migration velocity analysis [5, 6] and diffraction imaging [7, 8].Fomel [1] and Hubral et al. [9] point out that velocity continuation is a wave propagation process. Instead of wavefronts propagating as a function of time, images propagate as a function of migration velocity. Recent work has extended the concept to heterogeneous and anisotropic velocity models [10–15]. To account for anisotropy, the seismic velocity model must become multiparameter. Consequentially, velocity continuation generalizes to a process of implementing image transformations caused by changes in multiple parameters rather than the single isotropic velocity alone.Accounting for azimuthal variations in seismic velocity results in better event focusing and improved imaging in such media [16]. Azimuthal variation in velocity has been shown to be an indicator of preferentially aligned vertical fractures [17], lateral heterogeneity [18], regional stress [19], or a combination of these factors. However, velocity analysis is commonly first performed on prestack common midpoint (CMP) gathers, where the geologic cause of any observed azimuthal velocity variation is ambiguous. Without the help of additional diagnostic gathers such as hybrid or cross-spread gathers [20], or an interpretive comparison between picked root mean square (RMS) and interval velocities [21], the cause of azimuthal variations in velocity can be identified only after migration.Azimuthal seismic imaging commonly requires iterations between velocity analysis and imaging. Residual azimuthal variations in traveltime after migration can be measured by using migration binning schemes which preserve both offset and azimuth information [22, 23]. After the first pass of (isotropic) migration, azimuthal variations in velocity are detected from residual moveout, which then provides the velocity model for anisotropic migration. Iterative processing flows that use these strategies are popular not only because they are fairly efficient and intuitive, but also because they can be implemented with minimal modification to existing software. However, iterative imaging flows cannot guarantee convergence to the correct or optimal velocity model [24]. Velocity continuation has the underlying strategy of performing velocity analysis and imaging simultaneously and can thus be used to directly find an optimal velocity model without iteration. Sicking et al. [19] have demonstrated the success of a similar strategy of using imaging as a velocity analysis tool for 3D multiazimuth reflection seismic data. Azimuthal velocity continuation can provide a theoretical framework for this approach. With these benefits as motivation, we extend time-domain velocity continuation to 3D, accounting for the case of azimuthally variable migration velocity.


## 2. Theory

The theory of velocity continuation formulates the connection between the seismic velocity model and the seismic image as a wavefield evolution process. In doing so, the process can be implemented in the same variety of ways as seismic migration. Seismic migration in its many forms is commonly derived starting at the wave equation, which is approximated by its time and amplitude components by the eikonal and transport equations, and if necessary, a system of ray tracing equations. Velocity continuation is derived in the opposite order [2]. Starting with a geometrical description of the image, a corresponding kinematic equation for traveltime is derived to describe how the image moves according to changes in imaging parameters. Subsequently, the kinematic equation is used to derive a corresponding wave equation, which describes the dynamic behavior of the image as an evolution through imaging parameter coordinates. This section outlines the key steps of this derivation, starting with a traveltime equation that permits azimuthal variations in velocity.Grechka and Tsvankin [25] truncate a two-dimensional Taylor series expansion for a generally inhomogeneous anisotropic media to derive the “NMO ellipse” moveout equation. Geometrically, the NMO ellipse model still assumes that events have hyperbolic moveout with offset, but it allows the velocity to change with azimuth. We start here by using the same truncated 2D Taylor series expansion to describe an azimuth-dependent traveltime equation for the summation surface of zero-offset time migration,
(1)T2(x,y,W)=4(τ2+(x-y)TW(x-y)),
where τ is the one-way vertical traveltime after migration, x is the (x1,x2) surface position of the zero-offset receiver in survey coordinates, y is the surface position of the point source image, and superscript T denotes transpose. The three independent elements of the symmetric slowness matrix,
(2)W=(W11W12W12W22),
have units of slowness squared, and the eigenvalues and eigenvectors of W determine the symmetry axes of the effective anisotropic medium [25]. In most common geologic situations, the eigenvalues of  W are positive [26], and (1) describes an elliptical-hyperbolic traveltime surface in 3D—hyperbolic in cross-section view and elliptical in map view. The fast and slow moveout velocities are aligned with the major and minor axes of this ellipse. W11 and W22 are the squared moveout slownesses along their respective survey coordinates, x1 and x2. The third parameter, W12, arises from observing the ellipse in the x1-x2 survey coordinates, which are generally rotated relative to its major and minor axes.The three-parameter moveout model of (1) is analytically convenient and practical, but the parameters themselves are not intuitive to interpret in terms of more common geophysical or geological parameters. However, some simple geometric observations can help convert the three elements of W into more intuitive measurements. If the ellipse happens to be aligned with the survey coordinates, W12=0. Finding the rotation angle which properly diagonalizes W therefore allows one to predict the orientation of the symmetry axes. This amounts to an eigenvalue problem, where the fast and slow velocities can be found as the eigenvalues and eigenvectors of W. The eigenvalues, Wfast and Wslow, of the slowness matrix can be found following [25],
(3)Wslow,fast=12[W11+W22±(W11-W22)2+4W122].
Since the eigenvalues have units of slowness squared, the smaller eigenvalue is Wfast=1/vfast2. One can solve for the angle β between the acquisition coordinates and the symmetry axes by using the formula found by [25],
(4)β=tan-1[W22-W11+(W22-W11)2+4W1222W12].
The eigenvalues can then be used together with β to solve for the zero-offset migration slowness S as a function of source-receiver azimuth θ:
(5)S2(θ)=Wslowcos⁡2(θ-β)+Wfastsin2(θ-β).Equations (3)–(5) allow one to convert the mathematically convenient parameters of W to more intuitive parameters, such as the fastest and slowest propagation velocities (Vfast,Vslow), the azimuth of the slowest velocity (β), and the percent anisotropy (σ=100×(1-Vslow/Vfast)). Alternatively, W can be converted into other common geophysical parameterizations. For example, Grechka and Tsvankin [25] show that once the effective parameters W have been converted to slowness as a function of azimuth by (5), they can be expressed in terms of horizontal transverse isotropy parameters as
(6)S2(θ)=1VP021+2δ(v)sin2(θ)1+2δ(v),
where δ(v) is the Thomsen-style parameter [27], introduced by Tsvankin [28], and VP0 is the vertical P-wave velocity.Conventionally, one assumes that (1) characterizes a particular event defined in image coordinates (x,τ), but one can also describe how that event would transform given a change in the image parameters W. Regardless of the velocity model, the traveltime T must remain unchanged between different images. From this observation, we arrive at the following set of conditions:
(7)∇xT2=(∂T2∂x1∂T2∂x2)=8τ∇xτ+8W(x-y)=0,∇WT2=(∂T2∂W11∂T2∂W12∂T2∂W12∂T2∂W22)=8τ∇Wτ+4(x-y)(x-y)T=0.
Combining and reducing these conditions yields a system of equations that are defined only in the image parameter coordinates:
(8)2∂τ∂W11+τ(W22(∂τ/∂x1)-W12(∂τ/∂x2))2(W122-W11W22)2=0,2∂τ∂W22+τ(W12(∂τ/∂x1)-W11(∂τ/∂x2))2(W122-W11W22)2=0,(9)2∂τ∂W12-2τ(W12(∂τ/∂x1)-W11(∂τ/∂x2))(W22(∂τ/∂x1)-W12(∂τ/∂x2))(W122-W11W22)2=0.The system of kinematic equations describing azimuthally anisotropic velocity continuation is then found by combining (8)–(9). In a vector notation, this becomes
(10)∇Wτ+τ2W-1∇xτ(∇xτ)TW-1=0,
where ∇x and ∇W are in the form given by (7).The method of characteristics [29] provides a link between a kinematic equation (such as (10)) and its corresponding wave-type equation. Fomel [2] demonstrates specifically how the method can be used to derive a velocity continuation wave equation from its kinematic counterpart. By first setting the characteristic surface condition,
(11)ψ=t-τ(x,W)=0,
and replacing τ with ψ and t, we obtain an alternative form of (10):
(12)ψt∇Wψ+t2W-1∇xψ(∇xψ)TW-1=0.
Equation (11) guarantees that the wavefronts of time-domain image wavefield P exist only where the arrival time τ is equal to the recorded time t at a given location. Now take both ξi and ξj to represent each of t, W11, W12, W22, x1, and x2. According to the method of characteristics, if Λij is the coefficient in front of (∂ψ/∂ξi)(∂ψ/∂ξj) from kinematic equation (12), then the corresponding wave equation will have the same coefficients Λij in front of each ∂2P/∂ξi∂ξj derivative. The time derivative ψt is equal to 1 given (11), and is included in the first term of (12) to facilitate the use of the method of characteristics. Then, by introducing Pxx as the spatial Hessian matrix of the wavefield,
(13)Pxx=(∂2P∂x12∂2P∂x1∂x2∂2P∂x2∂x1∂2P∂x22),
we arrive at the azimuthally anisotropic post-stack velocity continuation wave equation,
(14)∇WPt=-t2W-1PxxW-1.
In the isotropic case, W is diagonal and W11=W22. Equation (14) then reduces to the isotropic velocity continuation equation first derived by Claerbout [30].


## 3. Method

Since velocity continuation is described by a wave equation, it can be implemented in analogous ways to seismic migration. Here, we demonstrate a spectral implementation of (14). By first stretching the time coordinate of an input image from t to t̃=t2/2, and then taking its 3D Fourier transform, (14) becomes the reduced partial differential equation:
(15)iΩ∇WP̂=12W-1kkTW-1P̂,
where Ω is the Fourier dual of t̃ and k is the wavenumber vector (Fourier dual of x). Equation (15) has the analytical solution,
(16)P̂(Ω,k1,k2,W)=P̂(Ω,k1,k2,W0)e-i/2ΩkT(W-1-W0-1)k,
which shows that continuation of an image from an arbitrary W0 to W can be achieved by multiplication with a shifting exponential in the Fourier domain. One can also directly migrate an unmigrated image by using the 2×2 matrix W0-1=0 for the initial velocity. In practice, the coordinate stretch from t to t̃ should be carefully applied as data will be compressed along the time axis for early samples.With a range of slowness matricesW, equation (16) can be used to quickly generate the corresponding range of anisotropically migrated images. When the correct velocity model is used, diffractions collapse to points, which we recognize as the image coming into focus. Although constant velocity models are used for each image, this type of spectral implementation can still be useful in the heterogeneous case, as different parts of the image will come into focus locally as the appropriate velocity is used [31, 32]. Once the range of images is generated, we search for the best focused image at each output location. We use the image attribute of kurtosis, defined as,
(17)ϕ(W)=∬P4(x,t,W)dxdt[∬P2(x,t,W)dxdt]2,
to quantify how well a location is focused in a particular image [8, 33]. Including integration limits specifies a window size for locally measuring kurtosis in the image. In application, either the integration limits control the size of a “sliding window,” or when viewing kurtosis as a local attribute, [34], they can be used to control the smoothness enforced by shaping regularization. In either case, the integration limits control a trade-off between the robustness of the focusing measurement and the resolution. From experience, typical limits for field data correspond to window sizes on the order of 101 samples in each dimension. It should be noted that the traveltime approximation of (1) loses accuracy in the presence of strong lateral heterogeneity, but is commonly used to estimate smooth effective parameter models. Following the maximum values through the resulting six-dimensional kurtosis hypercube, ϕ(t,x,W), and then slicing corresponding pieces from the output images volume, P(t,x,W), reveals an effective medium-based heterogeneous velocity model and a well focused image. This spectral implementation and slicing procedure is similar to searching through a set of constant-velocity f–k migrations and can be completed without any prior knowledge of the velocity model [32, 35].


## 4. Examples

Two simple synthetic examples are provided below to illustrate 3D velocity continuation over a range of velocity models. In the first example, we apply velocity continuation to a point diffractor. In the second example, we apply the method to a synthetic post-stack image of a set of faults. The second example illustrates fracture characterization through diffraction imaging as a potential application for 3D azimuthal velocity continuation. The data in both examples are modeled using (1), which geometrically approximates a diffraction surface as an elliptical-hyperbolic surface. Field data and more accurately modeled data will generally also exhibit nonhyperbolic moveout, for which (1) does not account. The physical validity and limitations of (1) are thoroughly discussed by Grechka and Tsvankin [25], but we focus here on how well diffractions can be collapsed, and how well the velocity parameters can be measured, if the data are ideally described by (1).Figure1(a) shows a single diffraction event, modeled using (1). The fastest direction of propagation is at β=105° counterclockwise from the x1 axis, with Vfast=3.50 km/s. The data in Figure 1(a) were modeled with σ=7% anisotropy, which may be quite high for most field cases, but it was chosen to allow the azimuthal variations in diffraction moveout to be visibly pronounced. As described above, we first stretch the time axis from t to t̃ and take the 3D Fourier transform of the data. Then we apply the phase-shift prescribed by (16) for a range of W. We found it more intuitive to specify the parameter ranges in terms of Vx1, β, and σ, and then convert them at each step into the three parameters of W for use in (16). The inverse of the in-line velocity squared 1/Vx12 is equivalent to W11, which, along with a given fast azimuth β and percent anisotropy σ, can be used to calculate W12 and W22 using (3)–(5). Last, we apply the 3D inverse Fourier transform and unstretch from t̃ to t to obtain the 6D image volume. Examples from the image volume using incorrect parameters are shown in Figures 1(b) and 1(c). The correct parameters are used in Figure 1(d), where the image is well focused.(a) A single azimuthally anisotropic diffraction. (b) The diffraction migrated by velocity continuation using correct parameters exceptσ=10, resulting in overmigration along x2. (c) Migration using the correct W11, but assuming isotropy. The result is now undermigrated along x2. (d) Migration using correct parameters. The image is well focused in both directions.
(a)
Diffraction(b)
Overmigrated(c)
Isotropic(d)
Correct parameterSince only a single diffraction is present in this example, we can measure kurtosis over a window spanning the entirety of each 3D image, reducing the kurtosis volume from 6D to 3D. Figure2 is a 2D slice of the kurtosis volume at the correct W11=1/Vx12 value of 0.0935 = 1/(3.27 km/s)2. The peak of the kurtosis map is near the correct values of σ=7 and β=105°. Once the peak of the kurtosis map is identified, one could refine the increments around the peak to yield more accurate estimates. The physical limitations of resolving azimuthal velocity parameters are discussed by Al-Dajani and Alkhalifah [37].Figure 2
Kurtosis values for the velocity continuation of the diffraction in Figure1(a). The map covers a range of anisotropy and angle values with an increment in β of 5° and an increment in σ of 0.5%. The correct values at 105° and 7% anisotropy (indicated by crosshairs) coincide with the peak of the kurtosis map.In practice, a conventional in-line 2D velocity analysis directly yieldsW11 from 1/Vx12, so Figure 2 could illustrate a realistic scenario for using 3D velocity continuation to improve upon a previous isotropic velocity model. In such a case, one would use previous Vx1 picks to hold W11 constant and then effectively test a variety of W12 and W22 values. Since W11 and W22 are measured with respect to the survey coordinates, either (or both) can be measured independently via a single-azimuth semblance scan, along x1 or x2, respectively. The best isotropic velocity based on a fully multiazimuth semblance scan will generally not represent either W11 or W22, but it can help limit the range of test parameters. Note that our method does not require prior knowledge of the velocity model, but without prior knowledge, the kurtosis measure remains a 6D volume. Although more difficult to visualize, the 6D kurtosis volume is computationally just as easily scanned for optimal imaging parameters as the 2D map in Figure 2.In the next example, we illustrate the concept of applying 3D anisotropic velocity continuation to diffraction imaging and fracture characterization. Figure3(a) shows a 3D synthetic post-stack diffraction data set, equivalent to the ideal separation of diffractions from specular reflections in post-stack data following Fomel et al. [8]. A fault map from Hargrove [36] (shown in Figure 3(a)) was digitized and used to create a 3D fracture model. Each fault location was used to generate a point diffraction in a homogeneous anisotropic background via (1). A timeslice of the modeled diffraction data is shown in Figure 3(b). The faults in the model typically have a strike of 112∘, and in cases where faults and nearby fractures (which more likely influence the seismic velocity) are similarly aligned, the fast direction of seismic wave propagation tends to align with their strike. By assuming a typical tight sandstone velocity of Vfast=4.0 km/s with 3% anisotropy, we choose the modeling W to be comprised of W11=0.0659, W22=0.0631, and W12=0.0014 (all in s2/km2). This results in a fast velocity direction along the strike of the faults. In Figure 3(d), we see that 3D velocity continuation using the correct parameters (again found by maximum kurtosis) allows the faults to be clearly imaged. If an intermediate isotropic velocity model is used, as in Figure 3(c), the diffractions are still imaged, but they are not as well focused compared to the anisotropically migrated diffractions in Figure 3(d). Conventionally, diffraction arrivals such as those in Figure 3(a) may be viewed as noise, but by separating them and treating them as signal, we can see here that imaging of steep and detailed features while simultaneously extracting anisotropy information may be possible.(a) Fault map from Northwest Scotland [36] used to model diffraction data. (b) Synthetic post-stack diffraction data modeled using (1) and a 3D model based on the fault map in (a). (c) Diffractions from (b) migrated using an isotropic velocity model. (d) Diffractions from (b) migrated by anisotropic 3D velocity continuation.
(a)
Fault map(b)
Fault diffractions(c)
Isotropic velocity continution(d)
Anisotropic continuation


## 5. Discussion and Conclusions

By extending time-domain velocity continuation to the azimuthally anisotropic 3D case, we have combined the concepts of azimuthal imaging and diffraction imaging. We assume a three-parameter migration slowness model that allows velocity to vary elliptically with azimuth. We have provided simple examples to illustrate the potential application of our method to fracture characterization through diffraction imaging. By treating diffractions as signal, our method performs azimuthal analysis on post-stack data, without the requirement for common-offset-vector or offset-vector-tile binning schemes. This is possible because, unlike reflections, diffractions can preserve azimuthal information post-stack. Post-stack data volumes have obvious advantages over prestack or vector-binned data for analysis, including smaller memory size and improved signal-to-noise ratio.Allowing azimuthal variation in the migration velocity will result in improved imaging, which is clearly a benefit of 3D velocity continuation. However, the potential for fracture characterization may be even more useful. Our method has many of the same ideas as the azimuthal imaging and fracture characterization flow proposed by Sicking et al. [19] for reflection data. Under the velocity continuation framework, we can extend the azimuthal imaging idea to 3D diffraction imaging. Since diffraction-generating fractures and faults are often nearly vertical and preferentially aligned, they can be associated with azimuthal anisotropy. Fomel et al. [8] demonstrate that it is possible to separate diffractions from specular reflections and then image their associated discontinuities through the use of velocity continuation. Their method operates on post-stack data, as they show that diffractions are highly sensitive to migration velocity, even in the zero-offset case. Al-Dajani and Fomel [38] have successfully demonstrated zero-offset diffraction image focusing as a fracture detection attribute on azimuth-sectored 3D field data. Our proposed method uses multiazimuth image focusing primarily as a velocity analysis criterion, but kurtosis could also be used as an image attribute. In cases where subsurface fractures cause azimuthal anisotropy, kurtosis as an attribute may be indicative of fracture intensity [38]. By applying velocity continuation to 3D diffraction imaging, one may be able to estimate both the orientation and the intensity of fractures from the resulting anisotropic velocity model and maximum kurtosis volumes, respectively. This information can be useful in reservoir development, as it can provide insight to subsurface fluid flow behavior.Although the spectral implementation of our method allows a range of possible images to be computed efficiently, it demands large amounts of memory to store a suite of images as well as the kurtosis volume. Modern computational hardware makes our approach feasible as is, especially for target-oriented imaging and analysis strategies. Future studies may lead to better optimization-based approaches to finding local kurtosis maxima, in which case, our method could be practical for dense parameter estimation throughout full 3D volumes.The underlying strategy of velocity continuation is to simultaneously estimate the velocity model as the data are imaged. This is beneficial in the case of azimuthal anisotropy discussed here, as the ambiguity between structural heterogeneity and anisotropy is handled without the need for iteration. Other multiparameter seismic imaging problems, such as converted-wave imaging, which are also conventionally handled by iterative flows, could also benefit from prestack versions of the 3D velocity continuation strategy.



---

*Source: 484653-2011-05-30.xml*