Quantifying and correcting residual azimuthal anisotropic moveout in image gathers using dynamic time warping

Introduction

The upper crust of the Earth is a complex and heterogeneous media containing different rock strata with constituent bedding planes, fracture networks, faults, and other features below seismic resolution that may cause the measurement of elastic properties to change with orientation. In such media, seismic waves propagate at different velocities depending on their direction of travel (Crampin, 1985,1981; Thomsen, 1988; Crampin, 1984a). This phenomena is referred to as seismic anisotropy, and there is a rich tradition of geophysical literature and research focused on its modeling and processing (Helbig and Thomsen, 2005; Grechka, 2009; Helbig, 1994; Thomsen, 2002; Tsvankin et al., 2010; Tsvankin, 2012; Crampin, 1984b).

Anisotropy type is determined by the symmetries present for wave propagation as a function of orientation. Transverse isotropy (TI) models media possessing a single axis of rotational symmetry (Thomsen, 1988; Crampin, 1986). Vertically transverse isotropy (VTI) refers to situations where horizontal wave propagation has the same velocity regardless of azimuthal orientation. In these circumstances velocity varies with inclination relative to bedding. This case is found to effectively model shales. Allowing the symmetry axis to tilt with dipping beds gives rise to tilted transverse isotropy (TTI). Horizontally transverse isotropy (HTI) occurs when the symmetry axis tilts fully to the horizontal. This type of anisotropy effectively models rocks with vertically aligned fracture networks and may be used to predict the fracture network orientation (Tod et al., 2007; Corrigan et al., 1996). Situations where velocity depends on both the inclination and azimuth of propagation lead to orthorhombic anisotropy and lower symmetry anisotropy systems. In the case of orthorhombic anisotropy the medium possesses three planes of symmetry rather than symmetry axes (Thomsen, 1988; Crampin, 1986; Tsvankin, 1997).

Orthorhombic anisotropy may be completely described by nine parameters, inverting for all of which becomes a computationally expensive exercise. This motivates the approximation of wave propagation with simpler anisotropy models possessing fewer parameters (Grechka et al., 2005), such as TI models which may be fully described with five (Tsvankin, 2012,1997). Stronger assumptions, namely that only vertical and horizontal velocities differ, allow approximation of p-wave anisotropy with just a single anisotropy parameter (Alkhalifah and Tsvankin, 1995). An approximation for HTI media wave propagation provided by Grechka and Tsvankin (1998) shows that the variation in velocity as a function of azimuth is elliptical. This type of HTI anisotropy is known as elliptical HTI anisotropy, and occurs when the phase slowness and group velocity surfaces are ellipsoidal. The kinematics of elliptical HTI anisotropy may be described by two anisotropy parameters plus velocity for each subsurface position (Abedi et al., 2019).

Seismic processing techniques incorporating anisotropy, which thus involve inverting for anisotropic parameters, are more computationally expensive than those that do not (Alkhalifah et al., 1996). Techniques that allow for more complex anisotropy, featuring additional parameters and thus more degrees of freedom, are still more expensive. Nonetheless, failing to fully account for anisotropy leads to seismic images that are less focused and accurate (Helbig and Thomsen, 2005). This is because many traces are migrated with what amounts to an incorrect velocity (Thomsen, 2001; Alkhalifah and Larner, 1994). The sensitivity of depth domain imaging techniques to velocity perturbation further accentuates these effects in depth images (Tsvankin et al., 2010). Seismic processing workflows must therefore balance the demands of accuracy and efficiency when treating anisotropy, and seismic processing practitioners are motivated to seek approximations that are able to account for anisotropy in a more efficient way (Helbig and Thomsen, 2005; Tsvankin et al., 2010). This desire prompted Burnett and Fomel (2009) to formulate a velocity - independent method for correcting azimuthal velocity variations using local travel time slopes in common midpoint (CMP) gathers. We propose to use dynamic time warping to inexpensively perform the correction.

The digital signal processing technique of dynamic time warping (DTW), developed by Sakoe and Chiba (1978) and applied to seismic imaging problems by Hale (2013), determines the set of integer shifts, $s [i ]$ , for a signal sample index $i$ , that most closely align a matching signal, $g [ i ]$ , to a reference signal, $f [ i ]$ , such that $f [ i ] \approx g [ i+s [ i ] ]$ . Constraints are applied to the process by declaring the maximum possible shift and the maximum strain, or how quickly shifts are permitted to change with respect to index $i$ . Shifts $s [i ]$ are determined by selecting a set of integer shifts obeying the imposed constraints that minimize the accumulated mismatch between the reference and matching signals over their entirety. Because only integer shifts are considered, calculation is a relatively rapid process.

When used in conjunction with a seismic migration method accounting for a VTI media, DTW enables us to correct for the residual azimuthal anisotropic moveout resulting from the un-accounted for elliptical HTI anisotropy in a computationally efficient manner. The algorithm is unconcerned with the physics of wave propagation and treats this moveout correction as a less expensive integer-shift data matching problem. If we further assume that the HTI fast axis is aligned with fracture networks present in the subsurface, determining the principal axes of the anisotropic ellipse with respect to wave propagation azimuth will provide us with a sort of average fracture network orientation over the whole ray path. Similarly, the elongation of this ellipse provides a notion of how anisotropic the material is over the whole ray path.

In the following sections this paper details how DTW may be used to create a method that compensates for residual elliptic HTI anisotropy in image gathers, resulting in enhanced seismic images and providing a measure of the orientation and relative intensity of that anisotropy, as well as the assumptions made in developing the method. The ability of the proposed technique to successfully recover high frequency signal and the principal anisotropic axis orientation is demonstrated on a synthetic gather featuring residual moveout caused by artificial elliptical HTI anisotropy. The method is then applied to a 3D field data set, generating sharper, more coherent images, as well as plausible information about HTI anisotropy. Finally, the strengths and limitations of the proposed approach are discussed, as well as promising avenues for future research directions.

Quantifying and correcting residual azimuthal anisotropic moveout in image gathers using dynamic time warping