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.
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, where is the one-way vertical traveltime after migration, is the surface position of the zero-offset receiver in survey coordinates, is the surface position of the point source image, and superscript denotes transpose. The three independent elements of the symmetric slowness matrix, have units of slowness squared, and the eigenvalues and eigenvectors of determine the symmetry axes of the effective anisotropic medium [25]. In most common geologic situations, the eigenvalues of 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. and are the squared moveout slownesses along their respective survey coordinates, and . The third parameter, , arises from observing the ellipse in the - 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 into more intuitive measurements. If the ellipse happens to be aligned with the survey coordinates, . Finding the rotation angle which properly diagonalizes 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 . The eigenvalues, and , of the slowness matrix can be found following [25], Since the eigenvalues have units of slowness squared, the smaller eigenvalue is . One can solve for the angle between the acquisition coordinates and the symmetry axes by using the formula found by [25], The eigenvalues can then be used together with to solve for the zero-offset migration slowness as a function of source-receiver azimuth :
Equations (3)–(5) allow one to convert the mathematically convenient parameters of to more intuitive parameters, such as the fastest and slowest propagation velocities (,), the azimuth of the slowest velocity (), and the percent anisotropy (). Alternatively, can be converted into other common geophysical parameterizations. For example, Grechka and Tsvankin [25] show that once the effective parameters have been converted to slowness as a function of azimuth by (5), they can be expressed in terms of horizontal transverse isotropy parameters as where is the Thomsen-style parameter [27], introduced by Tsvankin [28], and is the vertical P-wave velocity.
Conventionally, one assumes that (1) characterizes a particular event defined in image coordinates (), but one can also describe how that event would transform given a change in the image parameters . Regardless of the velocity model, the traveltime must remain unchanged between different images. From this observation, we arrive at the following set of conditions: Combining and reducing these conditions yields a system of equations that are defined only in the image parameter coordinates: The system of kinematic equations describing azimuthally anisotropic velocity continuation is then found by combining (8)–(9). In a vector notation, this becomes where and 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, and replacing with and , we obtain an alternative form of (10): Equation (11) guarantees that the wavefronts of time-domain image wavefield exist only where the arrival time is equal to the recorded time at a given location. Now take both and to represent each of , , , , , and . According to the method of characteristics, if is the coefficient in front of from kinematic equation (12), then the corresponding wave equation will have the same coefficients in front of each derivative. The time derivative 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 as the spatial Hessian matrix of the wavefield, we arrive at the azimuthally anisotropic post-stack velocity continuation wave equation, In the isotropic case, is diagonal and . 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 to , and then taking its 3D Fourier transform, (14) becomes the reduced partial differential equation: where is the Fourier dual of and is the wavenumber vector (Fourier dual of ). Equation (15) has the analytical solution, which shows that continuation of an image from an arbitrary to can be achieved by multiplication with a shifting exponential in the Fourier domain. One can also directly migrate an unmigrated image by using the matrix for the initial velocity. In practice, the coordinate stretch from to should be carefully applied as data will be compressed along the time axis for early samples.
With a range of slowness matrices , 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, 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 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, , and then slicing corresponding pieces from the output images volume, , 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 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).
Figure 1(a) shows a single diffraction event, modeled using (1). The fastest direction of propagation is at counterclockwise from the axis, with km/s. The data in Figure 1(a) were modeled with 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 to and take the 3D Fourier transform of the data. Then we apply the phase-shift prescribed by (16) for a range of . We found it more intuitive to specify the parameter ranges in terms of , , and , and then convert them at each step into the three parameters of for use in (16). The inverse of the in-line velocity squared is equivalent to , which, along with a given fast azimuth and percent anisotropy , can be used to calculate and using (3)–(5). Last, we apply the 3D inverse Fourier transform and unstretch from to 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) Diffraction
(b) Overmigrated
(c) Isotropic
(d) Correct parameter
Since 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. Figure 2 is a 2D slice of the kurtosis volume at the correct value of 0.0935 = 1/(3.27 km/s)2. The peak of the kurtosis map is near the correct values of and . 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].
In practice, a conventional in-line 2D velocity analysis directly yields from , 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 picks to hold constant and then effectively test a variety of and values. Since and are measured with respect to the survey coordinates, either (or both) can be measured independently via a single-azimuth semblance scan, along or , respectively. The best isotropic velocity based on a fully multiazimuth semblance scan will generally not represent either or , 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. Figure 3(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 km/s with 3% anisotropy, we choose the modeling to be comprised of , , and (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
(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.
Acknowledgments
The authors would like to thank Vladimir Bashkardin, Hesam Kazemeini, Siwei Li, Yang Liu, and Xiaolei Song for helpful suggestions and Peter Eichhubl and Pete Hargrove for providing the example fracture map. The authors would also like to thank the three anonymous reviewers for comments helpful in improving this paper. Funding for this project is provided by RPSEA through the “Ultra-Deepwater and Unconventional Natural Gas and Other Petroleum Resource” program authorized by the US Energy Policy Act of 2005. RPSEA (http://www.rpsea.org/) is a nonprofit corporation whose mission is to provide a stewardship role in ensuring the focused research, development and deployment of safe and environmentally responsible technology that can effectively deliver hydrocarbons from domestic resources to the citizens of the United States. RPSEA, operating as a consortium of premier US energy research universities, industry, and independent research organizations, manages the program under a contract with the US Department of Energy’s National Energy Technology Laboratory.