3D velocity-independent elliptically-anisotropic moveout correction

Theory

Following (Grechka and Tsvankin, 1998), the elliptical NMO equation in 3D can be written with the help of a truncated Taylor series expansion as,

$$t_0=\sqrt{t^2-\left(W_xx^2+W_yy^2+2W_{xy}xy\right)},$$ (1)

where $t$ is the event arrival time, $t_0$ is the moveout corrected time, and $x$ and $y$ are the components of full offset in the x and y survey directions, respectively. $W_x$ and $W_y$ are the conventionally-measured moveout slownesses squared values (along the same survey coordinates). Equation 1 describes NMO with elliptical velocity, where the third parameter, $W_{xy}$ , arises from observing the ellipse from rotated coordinates. In practice, one can perform elliptically anisotropic NMO using conventional velocity picking in the inline and crossline directions, but the principal directions of moveout must also be estimated. Rewriting equation 1 as a matrix-vector multiplication between the offset vector and the slowness matrix $\tensor{W}$ allows one to solve for the angle between the acquisition coordinates and the medium symmetry axes, denoted as $\alpha$ here. This can be done either by finding the eigenvectors of the system (Grechka and Tsvankin, 1998), or by using geometric arguments and well-known relations between the formulas for a rotated ellipse and its unrotated equivalent (Weisstein, 2009). Using the latter approach gives an expression for $\alpha$ in terms of conventional slowness parameters,

HTINMO
Figure 1. Elliptic NMO velocity depends on the natural symmetry axes (a-b), not the acquisition coordinates (x-y). The angle $\alpha$ describes the counter-clockwise rotation between the two coordinate frames.

$$\alpha =\frac{1}{2}\tan ^{-1}\left(\frac{2W_{xy}}{W_x-W_y}\right).$$ (2)

In this expression, $\alpha$ is the angle from a survey axis measured counter-clockwise toward the nearest symmetry axis. If $W_x$ is equal to $W_y$ , then the arc-tangent argument goes to infinity, corresponding to $\alpha =45^\circ$ . Although equation 2 is a straightforward way of finding the coordinate rotation angle, finding the eigenvalues and eigenvectors allows one to resolve between the fast and slow principal moveout directions. The eigenvalues $\lambda _1$ and $\lambda _2$ of the slowness matrix,

$$\tensor{W}=\begin{pmatrix}W_x & W_{xy} \cr W_{xy} & W_y \end{pmatrix},$$ (3)

can be found following Grechka and Tsvankin (1998). Rewritten here in our notation, they demonstrate how the eigenvalues,

$$\lambda _{1,2}=\frac{1}{2}\left[ W_x+W_y\pm \sqrt{(W_x-W_y)^2+4W_{xy}^2}\right],$$ (4)

can be used together with $\alpha$ to solve for the NMO slowness $S$ as a function of source-receiver azimuth $\theta$ :

$$S_{nmo}^2(\theta )=\lambda _1\cos^2(\theta -\alpha )+\lambda _2\sin^2(\theta -\alpha).$$ (5)

Since equation 1 describes the predicted hyperbolic traveltime curve on a CMP gather for a variety of cases where either real or apparent elliptical anisotropy may be present, the local slope of an event on an inline or crossline CMP gather can be related to the conventional moveout slowness parameters by taking the derivative of 1 with respect to $x$ and $y$ . Ignoring higher order terms and assuming the parameters vary slowly along $x$ and $y$ , gives a first-order approximation of how the measured slopes relate to conventional moveout parameters:

$$p_x(t,x,y)=\frac{\partial t}{\partial x}=\frac{W_xx+W_{xy}y}{t},$$ (6)

$$p_y(t,x,y)=\frac{\partial t}{\partial y}=\frac{W_yy+W_{xy}x}{t}.$$ (7)

By substitution back into 1, we arrive at the velocity-independent expression for 3D elliptical moveout in terms of local slopes:

$$t_0=\sqrt{t^2-t \left(p_xx+p_yy\right)}.$$ (8)

Notice that only two parameters ($p_x$ and $p_y$ ) must be measured to completely predict the NMO corrected time. More importantly, these parameters can be measured automatically using a local slope estimation algorithm, such as plane-wave destruction (Fomel, 2002). Equation 8 is a 3D extension for the 2D equation from Ottolini (1983).

Automated processes allow one to save time spent on a project, but it may seem that the insight and information gained during a more interactive conventional processing flow would be lost. A significant part of production velocity analysis involves picking or examining the velocity model directly, which provides an early and intuitive link between the seismic data and the subsurface geology. The velocity model and anisotropy information are themselves invaluable sources of geologic information. They also control the positioning of events in the final image, so an ability to extract these parameters is desirable.

The relation between local slopes and moveout velocity has been documented for the 2D case (Wolf et al., 2004; Ottolini, 1983; Fomel, 2007). In the 3D case where apparent azimuthal anisotropy is present, at least three conventional slowness or velocity-like values are needed to characterize moveout ($W_x$ , $W_y$ , and $W_{xy}$ from equation 1). Although equation 8 suggests they are not necessary for moveout in terms of local slopes, these values can be used to characterize anisotropy, and may also be useful for other subsequent processing. Simply rearranging equations 6 and 7 gives expressions for $W_x$ and $W_y$ :

$$W_x=\frac{tp_x-W_{xy}y}{x},$$ (9)

and,

$$W_y=\frac{tp_y-W_{xy}x}{y}.$$ (10)

Both of these parameters require an estimate of $W_{xy}$ . A first-order approximation of $W_{xy}$ can be found by differentiating equation 6 with respect to $y$ or equation 7 with respect to $x$ :

$$W_{xy}=t\frac{\partial p_x}{\partial y}+p_xp_y=t\frac{\partial p_y}{\partial x}+p_xp_y=tp_{xy}+p_xp_y.$$ (11)

Since slopes are measured as a local attribute, the inline and crossline local slopes comprise data volumes with the same dimensions and coordinates as the input CMP. Applying a 1D derivative filter to these volumes allows one to obtain either mixed-derivative in equation 11, and solve for the apparent anisotropy angle $\alpha$ , using equation 1. This angle can also be expressed in terms of local slopes. Combining equations 9, 10, and 11 yields,

$$W_x-W_y=\frac{t}{xy}[yp_x-xp_y+(p_{xy}+(x^2-y^2)\frac{p_xp_y}{t})].$$ (12)

Now everything needed to express $\alpha$ independently of velocity is found in equations 11 and 12. Combining them with equation 2 gives,

$$\alpha(t,x,y)=\frac{1}{2}\tan ^{-1}\left(\frac{2xy(tp_{xy}+p_xp_y)}{t(yp_x-xp_y)+(tp_{xy}+p_xp_y)(x^2-y^2)}\right).$$ (13)

Implementing equation 13 creates an attribute for each input data sample describing the counter-clockwise azimuthal angle between the symmetry coordinates and the acquisition coordinates. Applying NMO to this attribute volume yields $\alpha (t_0,x,y)$ , which should then theoretically be constant at each time-slice if the moveout were exactly described by equation 1.

Finding local estimates of slowness and anisotropy parameters using equations 9-13 remains at this point only an interesting theoretical idea. A more robust and practical approach to extracting velocity and anisotropy parameters is to exploit the shear number of volumetric slope measurements made to perform the velocity-independent NMO correction. For a given CMP with dimensions ( $n_x\times n_y\times n_t$ ), the NMO correction applies a shift of time-squared,

$$t_\Delta (t_0,x,y) =t^2(x,y)-t_0^2(x,y),$$ (14)

which can be automatically computed for every output coordinate using equation (8) and stored as another volume of the same dimensions. Once NMO is applied, the time axis of the CMP gather represents $t_0$ , so the slowness matrix $\tensor{W}$ and $\alpha$ should each be constant for a given time value. Each time-slice from either the data or one of the attribute volumes can be viewed as an ( $n_x \times n_y$ ) matrix, which can be re-indexed into a vector of length ( $n_x \times n_y$ ). If the x and y indexes from the time-slice are $i$ and $j$ respectively, then the value from position ($i,j$ ) in the matrix is mapped to the $k=i+jn_x$ position in the vector. Using this notation, a highly overdetermined problem follows from writing equation 1 as a matrix-vector multiplication:

$$\mathbf{t}_\Delta=\tensor{X}\mathbf{w}.$$ (15)

where the $k^{th}$ element of $\mathbf{t}_\Delta$ is,

$$t_{\Delta k}=t_\Delta (t_0,x_i,y_j),$$ (16)

the $k^{th}$ row of $\tensor{X}$ is given by the vector,

$$\mathbf{x}_{k}=\begin{pmatrix}x^2_i & y^2_j & 2x_iy_j \end{pmatrix},$$ (17)

and

$$\mathbf{w}=\begin{pmatrix}W_x \cr W_y \cr W_{xy} \end{pmatrix}.$$ (18)

Linear system 15 has ( $n_x \times n_y$ ) equations with only three unknowns. By solving 15 for each time-slice in the output CMP, we construct the slowness matrix $\tensor{W}(t_0)$ , and use it with equations 13, 4, and 5 to extract the coordinate rotation angle $\alpha (t_0)$ and the NMO slowness as a function of azimuth, $S(t_0,\theta )$ .

3D velocity-independent elliptically-anisotropic moveout correction