NMO-velocity surfaces and Dix-type formulae in anisotropic heterogeneous media

TL;DR

This paper introduces NMO-velocity surfaces in 3D anisotropic media, enabling Dix-type averaging and differentiation, which improves modeling and tomography of complex layered subsurface structures.

Contribution

It develops a novel formalism of NMO-velocity surfaces and associated algorithms for anisotropic heterogeneous media, extending traditional methods to 3D anisotropic cases.

Findings

NMO-velocity surfaces can be visualized as 3D plots of NMO velocity in all directions.

Dix-type averaging of NMO-velocity surface cross-sections models NMO ellipse and moveout.

The method is efficient for layered anisotropic media and suitable for anisotropic stacking-velocity tomography.

Abstract

Reflection moveout of pure modes recorded on conventional-length spreads is described by a normal-moveout (NMO) velocity that depends on the orientation of the common-midpoint (CMP) line. Here, we introduce the concept of NMO-velocity surfaces, obtained by plotting the NMO velocity as the radius-vector along all possible directions in 3-D space, and use it to develop Dix-type averaging and differentiation algorithms in anisotropic heterogeneous media. The intersection of the NMO-velocity surface with the horizontal plane represents the NMO ellipse that can be estimated from wide-azimuth reflection data. We demonstrate that the NMO ellipse and conventional-spread moveout as a whole can be modeled by Dix-type averaging of specifically oriented cross-sections of the NMO-velocity surfaces along the zero-offset reflection raypath. This formalism is particularly simple to implement for a stack of homogeneous anisotropic layers separated by plane dipping boundaries. Since our method involves computing just a single (zero-offset) ray for a given reflection event, it can be efficiently used in anisotropic stacking-velocity tomography.

Equations

$$V_{\mathrm{nmo}}^{-2}(\mathcal{L})=\mathcal{L}\mathbf{U}{\mathcal{L}}^{\mathbf{T}},$$

$$U_{km}={\tau }_0\frac{{\partial }^2\tau (\mathbf{x})}{\partial x_k\partial x_m}\equiv {\tau }_0\frac{\partial p_k(\mathbf{x})}{\partial x_m},(k,m=1,2,3).$$

$$V_{\mathrm{nmo}}^{-2}(\alpha )=U_{11}{\cos }^2\alpha +2U_{12}\sin \alpha \cos \alpha +U_{22}{\sin }^2\alpha .$$

$$W_{ij}={\tau }_0\frac{\partial p_i(\mathbf{x})}{\partial x_j}=U_{ij},(i,j=1,2).$$

$${\bf U}=\left(\begin{array}[]{ccc}W_{11}&W_{12}&U_{13}\\ W_{12}&W_{22}&U_{23}\\ U_{13}&U_{23}&U_{33}\end{array}\right).$$

$$\tau (\mathbf{x})=\frac{1}{\nu }{\cosh }^{-1}\left[1+\frac{{\nu }^2(x_1^2+x_2^2+x_3^2)}{2V_0(V_0+\nu x_3)}\right].$$

$$F(\mathbf{p},\mathbf{x})=0,$$

$${\bf U}=\left(\begin{array}[]{ccc}W_{11}&W_{12}&q_{,1}W_{11}+q_{,2}W_{12}-\tau_{0}\,F_{x_{1}}/F_{q}\\ \bullet&W_{22}&q_{,1}W_{12}+q_{,2}W_{22}-\tau_{0}\,F_{x_{2}}/F_{q}\\ \bullet&\bullet&q_{,1}^{2}W_{11}+2q_{,1}q_{,2}W_{12}+q_{,2}^{2}W_{22}-\tau_{0}\,(q_{,1}F_{x_{1}}+q_{,2}F_{x_{2}}+F_{x_{3}})/F_{q}\end{array}\right)\,.$$

$${\bf U}^{\rm hom}=\left(\begin{array}[]{ccc}W_{11}&W_{12}&q_{,1}W_{11}+q_{,2}W_{12}\\ \bullet&W_{22}&q_{,1}W_{12}+q_{,2}W_{22}\\ \bullet&\bullet&q_{,1}^{2}W_{11}+2q_{,1}q_{,2}W_{12}+q_{,2}^{2}W_{22}\end{array}\right).$$

$$q_{,1}{U}_{k1}^{\mathrm{hom}}+q_{,2}{U}_{k2}^{\mathrm{hom}}=U_{k3}^{\mathrm{hom}},(k=1,2,3).$$

$$\det {\mathbf{U}}^{\mathrm{hom}}=0.$$

$$\left\{\begin{array}[]{l}W_{11}\,\dfrac{e_{1}}{e_{3}}+W_{12}\,\dfrac{e_{2}}{e_{3}}=-q_{,1}W_{11}-q_{,2}W_{12}\,,\\ \\ W_{12}\,\dfrac{e_{1}}{e_{3}}+W_{22}\,\dfrac{e_{2}}{e_{3}}=-q_{,1}W_{12}-q_{,2}W_{22}\,.\end{array}\right.$$

$$\frac{e_i}{e_3}=-q_{,i},(i=1,2).$$

$$\frac{g_i}{g_3}=-q_{,i},(i=1,2).$$

$${\bf W}={{p_{1}q_{,1}+p_{2}q_{,2}-q}\over{q_{,11}q_{,22}-q_{,12}^{2}}}\left(\begin{array}[]{r r}q_{,22}&-q_{,12}\\ -q_{,12}&q_{,11}\end{array}\right),$$

$$q_{,ij}=-\frac{F_{p_i{p}_j}+F_{p_{i}q}{q}_{,j}+F_{p_{j}q}{q}_{,i}+F_{qq}{q}_{,i}{q}_{,j}}{F_q},(i,j=1,2),$$

$$V_{\mathrm{nmo},P}=V_{P0}\sqrt{1+2\delta }$$

$$\eta \equiv \frac{ϵ-\delta }{1+2\delta },$$

$$U_{11}^{\mathrm{VTI}}\equiv W_{11}^{\mathrm{VTI}}=\frac{1}{V_{\mathrm{nmo},P}^2}-p_1^2\left[1+2\eta (4y^2-9y+6)\right],$$

$$U_{12}^{\mathrm{VTI}}\equiv W_{12}^{\mathrm{VTI}}=0,$$

$$U_{22}^{\mathrm{VTI}}\equiv W_{22}^{\mathrm{VTI}}=\frac{1}{V_{\mathrm{nmo},P}^2}-2\eta p_1^2(2-y),$$

$$U_{13}^{\mathrm{VTI}}=-\frac{p_1\sqrt{1-y}}{V_{\mathrm{nmo},P}}\left[1+\delta -\eta y\frac{8y^2-15y+8}{1-y}\right],$$

$$U_{23}^{\mathrm{VTI}}=0,$$

$$U_{33}^{\mathrm{VTI}}=p_1^2\left[1+2\delta -4\eta y(1-2y)\right],$$

$${\tau }_0{\mathbf{U}}^{-1}=\frac{\partial \mathbf{x}}{\partial \mathbf{p}}.$$

$$\tau_{0}\,\left[{\bf W}^{\cal P}\right]^{-1}=\frac{\partial{\bf x}^{\cal P}}{\partial{\bf p}}\,\Biggl{|}_{{\bf p}^{\cal P}}\,.$$

$$\tau_{0,\ell}\,\left[{\bf W}^{\cal P}_{\ell}\right]^{-1}=\frac{\partial{\bf x}_{\ell}^{\cal P}}{\partial{\bf p}}\,\Biggl{|}_{{\bf p}^{\cal P}}\,.$$

$${\left[{\mathbf{W}}^{\mathcal{P}}\right]}^{-1}\sum _{\ell }{\tau }_{0,\ell }=\sum _{\ell }{\tau }_{0,\ell }{\left[{\mathbf{W}}_{\ell }^{\mathcal{P}}\right]}^{-1}.$$

$$\frac{d\mathbf{x}}{d{\tau }_0}=\frac{\partial H}{\partial \mathbf{p}}$$

$$\frac{d\mathbf{p}}{d{\tau }_0}=-\frac{1}{2}\frac{\partial H}{\partial \mathbf{x}},$$

$${\mathbf{W}}^{\mathcal{P}({\tau }_0+\Delta {\tau }_0)}({\tau }_0+\Delta {\tau }_0)=\mathbf{U}({\tau }_0+\Delta {\tau }_0)\bigcap \mathcal{P}({\tau }_0+\Delta {\tau }_0).$$

$${\left[\mathbf{W}(\ell +1)\right]}^{-1}=\frac{{\tau }_0(\ell ){\left[\mathbf{W}(\ell )\right]}^{-1}+{\tau }_{0,\ell +1}{\left[{\mathbf{W}}_{\ell +1}\right]}^{-1}}{{\tau }_0(\ell +1)},$$

$$\mathbf{x}\equiv [x_1,x_2,x_3]=h\mathcal{L}\equiv h[{\mathcal{L}}_1,{\mathcal{L}}_2,{\mathcal{L}}_3].$$

Full text

The intersection of the NMO-velocity surface with the horizontal plane represents the NMO ellipse that can be estimated from wide-azimuth reflection data. We demonstrate that the NMO ellipse and conventional-spread moveout as a whole can be modeled by Dix-type averaging of specifically oriented cross-sections of the NMO-velocity surfaces along the zero-offset reflection raypath. This formalism is particularly simple to implement for a stack of homogeneous anisotropic layers separated by plane dipping boundaries. Since our method involves computing just a single (zero-offset) ray for a given reflection event, it can be efficiently used in anisotropic stacking-velocity tomography.

Application of the Dix-type averaging to layered transversely isotropic media with a vertical symmetry axis (VTI) shows that the presence of dipping interfaces above the reflector makes the $P$-wave NMO ellipse dependent on the vertical velocity and anisotropic coefficients $\epsilon$ and $\delta$. In contrast, $P$-wave moveout in VTI models with a horizontally layered overburden is fully controlled by the NMO velocity of horizontal events and the Alkhalifah-Tsvankin coefficient $\eta\approx\epsilon-\delta$. Hence, in some laterally heterogeneous, layered VTI models $P$-wave reflection data may provide enough information for anisotropic depth processing.

I Introduction

NMO velocity estimated from reflection traveltimes recorded in common-midpoint geometry provides valuable information about the velocity field and anisotropic parameters of the subsurfaceTsvankin (2001). Although the relationship between the measured moveout velocity and elastic parameters becomes rather complicated if the model is heterogeneous and anisotropic, the azimuthal dependence of NMO velocity has a simple explicit form. Grechka and TsvankinGrechka and Tsvankin (1998) examined pure-mode reflection traveltimes recorded at a fixed CMP location along different azimuths $\alpha$ in the horizontal plane and showed that the NMO velocity $V_{\rm nmo}(\alpha)$ typically varies as an *ellipse * in the horizontal plane. [$V_{\rm nmo}(\alpha)$ may have a different form only if CMP traveltime decreases (i.e., reverse moveout) with offset in one or more directions.] In the special case of a homogeneous isotropic layer, the NMO ellipse was first obtained by LevinLevin (1971).

The orientation and semi-axes of the NMO ellipse depend on the spatial derivatives of the slowness vector at the CMP location. The simplicity and generality of this result arise because NMO velocity governs the wavefront curvature at zero offsetShah (1973); therefore, its azimuthal variation has to be a quadratic function in the spatial coordinates. Grechka, Theophanis and TsvankinGrechka et al. (1999a) extended the equation of the NMO ellipse to mode-converted waves in horizontally layered anisotropic models with a horizontal symmetry plane in each layer.

The elliptical azimuthal dependence of the NMO-velocity function was used to develop efficient algorithms for azimuthal stacking-velocity analysis and moveout correction in wide-azimuth 3-D surveysCorrigan et al. (1996); Grechka and Tsvankin (1999a). Even more importantly, the equation of the NMO ellipse provides a foundation for moveout inversion in arbitrary anisotropic media. For models with a horizontally layered overburden above a dipping reflector, the NMO ellipse at the surface represents a Dix-type average of the interval NMO ellipsesGrechka et al. (1999b). This equation, generalizing the classical DixDix (1955) result and its extensions for isotropic mediaShah (1973); Hubral and Krey (1980), can be used to reconstruct interval NMO ellipses from surface data and then invert them for the anisotropic parameters. The parameter-estimation methodology based on this approach was successfully implemented for several common anisotropic models including orthorhombic mediaGrechka and Tsvankin (1999b) and transverse isotropy with a verticalGrechka and Tsvankin (1998), horizontalContreras et al. (1999) and tiltedGrechka and Tsvankin (2000) symmetry axis.

The papers on parameter estimation listed above, however, consider only laterally homogeneous models above the reflector. Although Grechka, Tsvankin and CohenGrechka et al. (1999b) outlined an approach for computing NMO ellipses in arbitrary heterogeneous media, their algorithm is purely numerical and is difficult to apply in interval parameter estimation. The correction of NMO ellipses for lateral velocity variation introduced by Grechka and TsvankinGrechka and Tsvankin (1999a) is restricted to horizontal layers with a horizontal symmetry plane.

Here, we relax the assumption of Grechka and TsvankinGrechka and Tsvankin (1998) that the CMP line belongs to the horizontal plane and examine the NMO velocity measured along an arbitrary direction ${\cal L}$ in 3-D space. (One can imagine, for instance, recording reflection arrivals along an oblique or vertical borehole.) If the vectors $V_{\rm nmo}({\cal L})$ are plotted from the CMP location, their ends form the NMO-velocity surface, while the NMO ellipse is the intersection of this surface with the horizontal plane.

Even though under normal circumstances we cannot count on measuring $V_{\rm nmo}$ along many different directions in space, this new theory naturally leads to a concise Dix-type representation of NMO ellipses in heterogeneous anisotropic media. As an example, we construct Dix-type formulae for a stack of homogeneous anisotropic layers separated by plane dipping interfaces. In the practically important case of non-elliptical VTI media, the influence of dipping interfaces above the reflector may make the $P$-wave NMO ellipses measured at the surface dependent on the interval vertical velocities and Thomsen’s coefficients $\epsilon$ and $\delta$, thus affording opportunities for reconstructing the model in depth using $P$-wave moveout data.

II NMO-velocity surfaces in heterogeneous anisotropic media

II.1 General formulation

We consider the NMO velocity of a pure-mode reflected wave that was recorded along an arbitrary oriented CMP line ${\cal L}$ in heterogeneous anisotropic media. It is assumed that the traveltime of the selected reflection event is uniquely defined for each (moderate compared to the reflector depth) source-receiver offset. If the traveltime becomes multi-valued, as in the vicinity of shear-wave cusps, the moveout function usually requires a more elaborate approximation than the hyperbolic equation parameterized by NMO velocity.

The exact function $V_{\rm nmo}({\cal L})$ is derived in Appendix A as [see equation (46)]

[TABLE]

where ${\cal L}=[{\cal L}_{1},\,{\cal L}_{2},\,{\cal L}_{3}]$ is a unit row vector, ${\cal L}^{\bf T}$ is a unit column vector, and ${\bf U}$ is a $3\times 3$ symmetric matrix with the elements

[TABLE]

Here $\tau_{0}$ is the one-way traveltime from the zero-offset reflection point to the CMP location, and $p_{k}({\bf x})$ are the components of the slowness vector ${\bf p({\bf x})}=[p_{1}({\bf x}),p_{2}({\bf x}),p_{3}({\bf x})]$ corresponding to rays excited at the zero-offset reflection point and recorded at location ${\bf x}$. The derivatives in equation (2) are evaluated at the common midpoint.

Azimuthally dependent NMO velocity in the horizontal plane (usually an ellipse) described by Grechka and TsvankinGrechka and Tsvankin (1998) can be viewed as the intersection of the NMO-velocity surface ${\bf U}$ with the horizontal plane. Substituting a horizontal unit vector ${\cal L}^{\rm hor}=[\cos\alpha,\,\sin\alpha,\,0]$ into the general expression (1) yields the NMO ellipse as a function of the azimuth $\alpha$,

[TABLE]

The $2\times 2$ matrix ${\bf W}$ introduced by Grechka and TsvankinGrechka and Tsvankin (1998) to define the NMO ellipse (3) coincides with the upper left submatrix of ${\bf U}$:

[TABLE]

II.2 Possible shapes of NMO-velocity surfaces

Equation (1) indicates that the function $V_{\rm nmo}({\cal L})$ defines a *centered quadratic * surface in 3-D space. The shape of this surface is determined by the eigenvalues of the matrix ${\bf U}$, which have to be real because ${\bf U}$ is real and symmetric. Using equation (4), ${\bf U}$ can be written in the form

[TABLE]

Equation (5) allows us to use the known properties of the matrix ${\bf W}$Grechka and Tsvankin (1998); Grechka et al. (1999b) to make several important observations about the matrix ${\bf U}$. First, if the NMO ellipse ${\bf W}$ has been found from moveout data, we need to determine only three quantities $U_{k3}$ to reconstruct the whole NMO-velocity surface ${\bf U}$. Below, we show that the elements $U_{k3}$ can be computed by differentiating the Christoffel equation at the CMP location. Second, since for homogeneous anisotropic media the elements $W_{ij}$ can be obtained in an explicit formGrechka et al. (1999b), the matrix ${\bf U}$ as a whole can be found explicitly as well. Third, the cross-section of the NMO-velocity surface ${\bf U}$ along the horizontal plane is elliptical because the matrix ${\bf W}$ typically represents an ellipse in the horizontal planeGrechka and Tsvankin (1998).

There are only three distinct types of quadratic surfaces which have elliptical cross-sections symmetric with respect to the CMP: ellipsoids, elliptical cylinders and one-sheeted hyperboloids, as shown in Figure 1. (Note that a hyperboloid and a cylinder may also have a non-elliptical intersection with the horizontal plane, as illustrated below.) Since the NMO-velocity surface is quadratic, it may have other shapes, such as those of a two-sheeted hyperboloid, an imaginary elliptical cylinder, a hyperbolic cylinder, etc. However, this may happen only if the matrix ${\bf U}$ has at least two non-positive eigenvalues, and reflection traveltime does not increase with offset in two or more directions in space. Although such cases are not prohibited by the theory, their occurrence is expected to be rare.

Note that if the NMO-velocity surface has the form of a cylinder, the NMO velocity along the axis of the cylinder is infinite, which implies that the traveltime in this direction does not change with offset. A numerical example below shows that both elliptical cylinders and one-sheeted hyperboloids can be encountered in realistic subsurface models; a more comprehensive discussion of cylindrical NMO-velocity surfaces is presented below.

II.3 Example of NMO-velocity surfaces

To investigate the shape of NMO-velocity surfaces for typical seismological models, we computed the matrix ${\bf U}$ for an isotropic medium with a constant vertical-velocity gradient [the velocity function is defined as $V(x_{3})=V_{0}+\nu x_{3}$]. For this model, the one-way traveltime $\tau({\bf x})$ from the origin of the coordinate system to point ${\bf x}=[x_{1},x_{2},x_{3}]$ can be found analyticallySlotnick (1959),

[TABLE]

Using equations (2) and (6), we derived explicit expressions for the matrix ${\bf U}$ in terms of $V_{0}$, $\nu$, and the depth and dip of the reflector.

Figure 1 displays the NMO surfaces, along with computed ray trajectories (circular arcs), for reflectors beneath this constant-gradient isotropic medium. Depending on reflector dip, the NMO-velocity surface can take any of the three shapes discussed above (an ellipsoid, a cylinder and a hyperboloid). Numerical tests for other isotropic models, where velocity monotonically increases with depth, indicate that NMO-velocity ellipsoids correspond to reflector dips below 90∘, cylinders to vertical reflectors, and hyperboloids to dips exceeding 90∘ (overhangs).

III Computation of NMO-velocity surfaces

III.1 Heterogeneous arbitrary anisotropic media

The matrix ${\bf U}$ contains six quantities [equation (5)] including three components of the matrix ${\bf W}$ responsible for the NMO ellipse. Computation of ${\bf W}$ was described by Grechka, Tsvankin and CohenGrechka et al. (1999b) and is further discussed below. The remaining elements $U_{3k}$ of the matrix ${\bf U}$ depend on the spatial derivatives of the vertical slowness component $p_{3}$, which can be obtained by combining $W_{ij}$ with the solution of the Christoffel equation at the common midpoint. Indeed, the slownesses at each spatial location ${\bf x}$ are related to each other by the Christoffel equation, which can be written as

[TABLE]

where ${\bf p}$ is the slowness vector at the spatial location ${\bf x}$ (e.g., Grechka, Tsvankin and Cohen, 1999). Equation (7) contains a separate contribution of the coordinates ${\bf x}$ because of the spatial dependence of the stiffness coefficients $c_{ij}$ in heterogeneous media. At a fixed location ${\bf x}$, equation (7) allows us to express the vertical slowness $p_{3}$ through the horizontal slownesses $p_{1}$ and $p_{2}$. Since the elements $U_{3k}$ depend on $p_{3}$, they can be found as functions of the derivatives $\partial p_{1}/\partial x_{k}$ and $\partial p_{2}/\partial x_{k}$ that determine the NMO ellipse ${\bf W}$. A complete derivation given in Appendix B leads to the following expression for the matrix ${\bf U}$ [equation (55)]:

[TABLE]

Here $q\equiv p_{3}$ is the vertical component of the slowness vector, $q_{,i}\equiv\partial q/\partial p_{i}=-F_{p_{i}}/F_{q}$, $F_{p_{k}}\equiv\partial F/\partial p_{k}$, $F_{q}\equiv\partial F/\partial q$ and $F_{x_{k}}\equiv\partial F/\partial x_{k}$. The matrix (8) is symmetric, so the bullets are used to denote the elements $U_{21}=U_{12}$, $U_{31}=U_{13}$ and $U_{32}=U_{23}$.

The derivatives $F_{x_{k}}$ with respect to the spatial coordinates in equation (8) depend on the medium properties (i.e., heterogeneity and anisotropy) near the CMP location. Therefore the NMO-velocity surface as a whole can be reconstructed from the NMO ellipse ${\bf W}({\bf x})$, if we know the slowness vector of the zero-offset ray and local values of the elastic constants $c_{ij}({\bf x})$ near the common midpoint.

III.2 Homogeneous media

Here we show that the NMO-velocity surface in homogeneous anisotropic media always represents a cylinder with the axis parallel to the zero-offset ray. The matrix ${\bf U}$ that describes this cylinder can be obtained in closed form using the Christoffel equation expressed in terms of the slowness components.

If the medium is homogeneous, the spatial derivatives $F_{x_{k}}$ vanish, and the matrix ${\bf U}$ from equation (8) simplifies to

[TABLE]

To find the shape of the corresponding NMO-velocity surface, note that the third column of the matrix ${\bf U}^{\rm hom}$ is a linear combination of the first two columns:

[TABLE]

As follows from equation (10),

[TABLE]

Since the first and the second columns are generally independent, the matrix ${\bf U}^{\rm hom}$ has one zero eigenvalue, so the surface defined by ${\bf U}^{\rm hom}$ has to be a cylinder. For models in which the matrix ${\bf W}$ describes an ellipse in the horizontal plane, the NMO-velocity surface is an elliptical cylinder. This conclusion is valid for any pure-mode reflections in homogeneous arbitrary anisotropic media.

The axis of the NMO-velocity cylinder is parallel to the eigenvector ${\bf e}=[e_{1},e_{2},e_{3}]$ corresponding to the zero eigenvalue of the matrix ${\bf U}^{\rm hom}$. Substituting the eigenvector ${\bf e}$ into the first two rows of the matrix (9) yields

[TABLE]

Therefore,

[TABLE]

The meaning of equation (13) can be explained using the expression for the components of the group-velocity vector ${\bf g}$ obtained by Grechka, Tsvankin and CohenGrechka et al. (1999b) [their equation (B–3)],

[TABLE]

Comparison of equations (13) and (14) shows that ${\bf e}$ is parallel to the vector ${\bf g}$ at the CMP location. In other words, the axis of the cylinder in homogeneous media of any symmetry points in the direction of the zero-offset ray. According to the geometrical meaning of the NMO-velocity surface, this result implies that the NMO velocity on the CMP line parallel to the zero-offset ray is infinite. Indeed, if sources and receivers are placed on the (straight) zero-offset ray, the reflected rays travel along the acquisition line; consequently, the two-way reflection traveltime in CMP geometry has to be independent of offset (i.e., $V_{\rm nmo}$ goes to infinity).

The exact equation of the NMO ellipse ${\bf W}$ in anisotropic homogeneous media is given by Grechka, Tsvankin and CohenGrechka et al. (1999b) [equation (7)] in terms of the slowness components,

[TABLE]

where

[TABLE]

$F_{p_{i}p_{j}}\equiv{\partial^{2}F}/{\partial p_{i}}{\partial p_{j}}$, $F_{p_{i}q}\equiv{\partial^{2}F}/{\partial p_{i}}{\partial q}$, and $F_{qq}\equiv{\partial^{2}F}/{\partial q^{2}}$. Since all quantities in equation (15) can be obtained explicitly from the Christoffel equation, equations (9) and (15) indicate that the whole NMO-velocity cylinder can also be constructed *analytically * for a given slowness vector ${\bf p}=[p_{1},p_{2},q]$ of the zero-offset ray.

III.3 NMO cylinder of the $\bm{P}$-wave in a weakly anisotropic VTI layer

According to equation (9), if the NMO-velocity cylinder ${\bf U}^{\rm hom}$ has been reconstructed from seismic data, it should be possible to find the derivatives $q_{,i}$ in addition to the NMO ellipse ${\bf W}$. Since for some models $q_{,i}$ may depend on medium parameters not constrained by the NMO ellipse, the NMO-velocity surface may provide valuable information for anisotropic inversion. This point is illustrated here for the $P$-wave NMO-velocity cylinder ${\bf U}^{\rm VTI}$ from a plane dipping reflector beneath a homogeneous VTI layer.

To simplify the derivation of the matrix ${\bf U}^{\rm VTI}$, we assume that the anisotropy is weak, and the NMO velocity can be linearized in Thomsen’sThomsen (1986) anisotropic coefficients $\epsilon$ and $\delta$. The $P$-wave NMO ellipse [i.e., the elements $W_{ij}$ in equation (9) expressed through the horizontal slowness components $p_{1}$ and $p_{2}$] in VTI media is fully controlled by the NMO velocity from a horizontal reflector

[TABLE]

and the anellipticity coefficient $\eta$ defined as

[TABLE]

where $V_{P0}$ is the $P$-wave vertical velocityAlkhalifah and Tsvankin (1995); Grechka and Tsvankin (1998). Hence, it is instructive to express our results in terms of $V_{{\rm nmo},P}$, $\eta$, and one of the generic Thomsen parameters (e.g., $\delta$) instead of the more conventional parameter set [$V_{P0}$, $\epsilon$, $\delta$].

Selecting the coordinate frame in which the reflector normal lies in the vertical plane $[x_{1},x_{3}]$, so that the zero-offset slowness component $p_{2}=0$ (Figure 2a), we use equations (9) and (15) to obtain

[TABLE]

where $y\equiv p_{1}^{2}\,V_{{\rm nmo},P}^{2}$.

Equations (18) – (20) describe the NMO ellipse ${\bf W}^{\rm VTI}$ with axes pointing in the dip and strike directions of the reflector. The dip component of the $P$-wave NMO velocity is defined by $W^{\rm VTI}_{11}$ in equation (18)Alkhalifah and Tsvankin (1995), while equation (20) for $W^{\rm VTI}_{22}$ gives the strike componentGrechka and Tsvankin (1998). Clearly, the NMO ellipse ${\bf W}^{\rm VTI}$ as a whole is governed by $V_{{\rm nmo},P}$ and $\eta$, with no dependence on $\delta$; this conclusion holds for strong anisotropy as wellGrechka and Tsvankin (1998).

Equations (21) – (23), specifying the additional components of ${\bf U}^{\rm VTI}$ needed to build the NMO-velocity cylinder, indicate that the $V_{\rm nmo}$ in non-horizontal directions depends on all three parameters ($V_{{\rm nmo},P}$, $\eta$ and $\delta$). This result also follows from the equation of the NMO ellipse in TI media with a horizontal symmetry axis (HTI) given inContreras et al. (1999). Note that the vertical symmetry axis ${\bf a}$ in Figure 2a becomes horizontal after rotating the whole plot by 90∘ around the coordinate axis $x_{2}$. Hence, this rotation transforms the VTI model in Figure 2a into the HTI model in Figure 2b. The quantities $U^{\rm VTI}_{22}$ and $U^{\rm VTI}_{33}$ determine the NMO ellipses in both the vertical $[x_{2},x_{3}]$-plane of the original VTI model and the horizontal $[x_{2},\tilde{x}_{3}]$-plane of the new HTI model. The results of Contreras et al.Contreras et al. (1999) show that the HTI ellipse ${\bf W}^{\rm HTI}$ is a function of $\delta$ as well as of $V_{{\rm nmo},P}$ and $\eta$.

Although the discussion of NMO velocities measured outside the horizontal plane may seem purely academic (unless vertical or oblique boreholes are available), intersections of the NMO-velocity surface with non-horizontal planes play an important role in Dix-type averaging of NMO velocities in heterogeneous anisotropic media. As we demonstrate below, the information contained in the matrix elements $U_{k3}$ can be extracted from the surface NMO ellipse in the presence of lateral heterogeneity above the reflector (such as intermediate dipping interfaces). It is known, though, that this is not possible for media with elliptical anisotropyDellinger and Muir (1988).

IV Dix-type formulae in heterogeneous anisotropic media

The NMO-velocity surfaces ${\bf U}({\bf x})$ can be called “effective” because, just as for effective NMO velocities, they incorporate the influence of the medium properties along the whole ray path between the zero-offset reflection point and the CMP location. Here, we devise Dix-type formulae for building the effective NMO-velocity surfaces from interval (or local) surfaces in heterogeneous anisotropic media.

IV.1 General considerations

Let us assume that the projection of the slowness vector ${\bf p}$ onto a certain plane ${\cal P}$ is preserved along the segment $L$ of the zero-offset ray. That will be the case, for example, if the ray crosses homogeneous layers separated by plane parallel interfaces ${\cal P}$. For simplicity, suppose that this segment starts at the zero-offset reflection point. To find the NMO-velocity surface ${\bf U}$ at the end of the segment $L$, it is convenient to rewrite equation (2) in the vector form,

[TABLE]

Constructing the intersection of the surface ${\bf U}$ with the plane ${\cal P}$, we obtain the NMO ellipse ${\bf W}^{\cal P}={\bf U}\bigcap{\cal P}$, which satisfies

[TABLE]

The superscript “${\cal P}$” in ${\bf x}^{\cal P}$ and ${\bf p}^{\cal P}$ emphasizes that the derivatives of the ray coordinates ${\bf x}$ are taken in the plane ${\cal P}$ $({\bf x}^{\cal P}={\bf x}\in{\cal P})$ with respect to the projection of the slowness vector ${\bf p}$ onto this planeGrechka et al. (1999b).

Next, we divide the segment $L$ into a number of smaller intervals $\ell$ and define the interval zero-offset traveltimes $\tau_{0,\ell}$ and NMO ellipses ${\bf W}^{\cal P}_{\ell}$. The interval NMO ellipses correspond to non-existent reflectors orthogonal to the slowness vector of the zero-offset ray. Applying equation (25) to each interval yields

[TABLE]

Summing up equation (26) over the segment $L$ and taking into account that ${\bf x}^{\cal P}=\sum_{\ell}{\bf x}^{\cal P}_{\ell}$ and $\tau_{0}=\sum_{\ell}\tau_{0,\ell}$, we find

[TABLE]

Equation (27) is identical to the generalized Dix formula of Grechka, Tsvankin and CohenGrechka et al. (1999b) derived for horizontally layered media (i.e., for a horizontal plane ${\cal P}$) above a dipping reflector. Note, however, that the plane ${\cal P}$ in equation (27) can have an arbitrary orientation.

The derivation above can be repeated for the segment $L$ located anywhere on the zero-offset ray. To obtain the intersection of the NMO-velocity surface ${\bf U}$ with the plane ${\cal P}$ at the end of $L$ we have to compute $\tau_{0}\,\left[{\bf W}^{\cal P}\right]^{-1}$ at the beginning of the segment and add it to the right-hand side of equation (27). Therefore, if the projection of the slowness vector ${\bf p}^{\cal P}$ onto the plane ${\cal P}$ is preserved along the segment $L$, the contribution of this segment to the intersection of the effective NMO-velocity surface with ${\cal P}$ can be obtained using Dix-type averaging of the corresponding intersections ${\bf W}^{\cal P}_{\ell}$ of the interval NMO-velocity surfaces ${\bf U}_{\ell}$.

Below, we demonstrate how the effective NMO velocity in heterogeneous anisotropic media can be obtained by integrating the interval NMO ellipses along the zero-offset ray.

IV.2 Heterogeneous anisotropic media

Suppose the medium is heterogeneous, and all components of the slowness vector vary in some fashion along the zero-offset ray. The description of rays in heterogeneous anisotropic media is given by the following system of differential equationsČervený et al. (1977):

[TABLE]

and

[TABLE]

where $\tau_{0}$ is the traveltime along the ray, and $H\equiv H({\bf p},{\bf x})$ is the Hamiltonian of some particular form, which does not need to be specified here. Equation (29) indicates that the slowness ${\bf p}$ changes in the direction $d{\bf p}/d\tau_{0}$ as we move along the ray. Hence, the projection of the slowness vector onto the tangent plane ${\cal P}(\tau_{0})\,\bot\,d{\bf p}/d\tau_{0}$ (i.e., the plane ${\cal P}(\tau_{0})$ is orthogonal to the vector $d{\bf p}/d\tau_{0}$) is locally preserved (Figure 3). Therefore, equation (27) can be applied to the NMO ellipse ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0})={\bf U}(\tau_{0})\bigcap{\cal P}(\tau_{0})$ at the infinitesimal ray segment corresponding to the interval traveltime $\Delta\tau_{0}$ to produce the ellipse ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0}+\Delta\tau_{0})$.

To account for the fact that ${\cal P}(\tau_{0}+\Delta\tau_{0})$ generally differs from ${\cal P}(\tau_{0})$ (Figure 3), we reconstruct the whole NMO surface ${\bf U}(\tau_{0}+\Delta\tau_{0})$ [using the Christoffel equation] and find its intersection with the plane ${\cal P}(\tau_{0}+\Delta\tau_{0})$ (see Appendices C and D):

[TABLE]

The resulting NMO ellipse ${\bf W}^{{\cal P}(\tau_{0}+\Delta\tau_{0})}(\tau_{0}+\Delta\tau_{0})$ can be continued along the next time interval, starting at $\tau_{0}+\Delta\tau_{0}$. Using this formalism, the NMO surface can be built by integrating the local NMO ellipses while solving the ray-tracing equations in heterogeneous anisotropic media. A more detailed mathematical description of this procedure is given in Appendices C and D. On the whole, the results above show that it is possible to model NMO ellipses in heterogeneous anisotropic media by tracing a single (zero-offset) ray.

IV.3 Homogeneous layers separated by plane dipping interfaces

The theory of the Dix-type averaging of NMO-velocity surfaces yields relatively simple results for the practically important special case of piecewise homogeneous media composed of anisotropic layers (or blocks) separated by plane dipping interfaces. In such a medium, the projection of the zero-offset slowness vector onto each interface is preserved due to Snell’s law (i.e., ${\bf p}_{\ell}\times{\bf z}_{\ell}={\bf p}_{\ell+1}\times{\bf z}_{\ell}$ at the $\ell$th interface with the normal ${\bf z}_{\ell}$; see Figure 4). Therefore, the layer boundaries play the role of the planes ${\cal P}$ that determine the intersections ${\bf W}^{\cal P}$ to be averaged by the Dix-type equation. Note that, as shown above [equation (9)], NMO-velocity surfaces in piecewise homogeneous media always have a cylindrical shape.

Figure 4 schematically illustrates the 3-D process of constructing the NMO-velocity cylinders in layered media with dipping interfaces. Assuming that the slownesses ${\bf p}_{\ell}$ and traveltimes $\tau_{0,\ell}$ have already been obtained from ray tracing, the Dix-type averaging can be performed as follows:

Step 1.

Using equations (9) and (15), compute the NMO-velocity cylinder ${\bf U}_{1}$ in the layer $\ell=1$ immediately above the reflector; the slowness vector ${\bf p}_{1}$ is parallel to the reflector normal. The interval cylinder ${\bf U}_{1}$ is equal to the “effective” cylinder ${\bf U}(1)$ in the first layer. If the layer number $\ell>1$, the cylinder ${\bf U}(\ell)$ (dashed lines in Figure 4a) is obtained from the continuation procedure described here.

Step 2.

Apply equation (63) to determine the intersection ${\bf W}(\ell)$ (the magenta line in Figure 4a) of the cylinder ${\bf U}(\ell)$ with the $\ell$th interface that has the normal ${\bf z}_{\ell}$.

Step 3.

Compute the interval cylinder ${\bf U}_{\ell+1}$ (dashed lines in Figure 4b) using equations (9) and (15) for the slowness vector ${\bf p}_{\ell+1}$. Find the intersection ${\bf W}_{\ell+1}$ (the magenta line in Figure 4b) of the cylinder ${\bf U}_{\ell+1}$ with the $\ell$th interface.

Step 4.

Obtain the cross-section ${\bf W}(\ell+1)$ of the effective cylinder ${\bf U}(\ell+1)$ at the top of the $\ell+1$th layer (the magenta line in Figure 4c) from the Dix-type formula (27):

[TABLE]

where $\tau_{0}(\ell)=\sum_{\jmath=1}^{\ell}\tau_{0,\jmath}$.

Step 5.

Reconstruct the cylinder ${\bf U}(\ell+1)$ (dashed lines in Figure 4d) using equations (68) and (9).

Step 6.

Repeat Step 2 for the next ($\ell+1$)th layer.

The sequence described above makes it possible to compute NMO ellipses for layered media with plane dipping interfaces without multi-azimuth, multi-offset ray tracing. Our Dix-type averaging procedure was used by Grechka, Pech and TsvankinGrechka et al. (2000a, b) to devise an efficient algorithm for $P$-wave stacking-velocity tomography in piecewise homogeneous VTI media.

IV.4 Numerical example

To verify the accuracy of Dix-type averaging in layered media with dipping interfaces, we computed the NMO ellipses at the horizontal surface for a model composed of three transversely isotropic layers with a tilted symmetry axis (Table 1). Figure 5 displays the NMO ellipses for the reflection from the bottom of the model determined from the Dix-type averaging procedure (solid) and 3-D anisotropic ray tracing (dotted). The ellipses almost coincide, thus confirming that the Dix-type equations give an adequate description of reflection moveout on conventional-length spreads. The small difference of up to 1.6% between the theoretical and ray-traced ellipses in Figure 5 can be attributed to the influence of nonhyperbolic moveout, which is not taken into account by our NMO-velocity equations. However, the errors due to nonhyperbolic moveout are small (at least, for $P$-waves), when the maximum offset does not exceed roughly the distance between the CMP and the reflector. This conclusion holds for $P$-wave data in a wide variety of anisotropic models of different complexityTsvankin (2001); Tsvankin and Thomsen (1994); Grechka and Tsvankin (1999b).

V Discussion

V.1 General results

This work introduces the concept of NMO-velocity surfaces in anisotropic heterogeneous media and applies the new theory to devise Dix-type averaging procedures for effective NMO velocity. Because the NMO-velocity surface is quadratic, it depends on six generally independent elements of a symmetric matrix $U_{km}$, which include both effective quantities averaged between the reflector and CMP location and local quantities defined at the common midpoint. If the medium near the CMP is homogeneous, the NMO-velocity surface always represents a cylinder, irrespective of the complexity of the model as a whole. Other shapes that may be encountered even in isotropic media are an ellipsoid and a one-sheeted hyperboloid.

The surfaces ${\bf U}$ provide the most general description of conventional-spread normal moveout because they can be used to determine NMO velocity in any direction in 3D space. One important practical example is the NMO ellipse formed by NMO velocities plotted in all possible azimuths within a certain plane. Since the NMO ellipse can be viewed as a cross-section of the NMO-velocity surface ${\bf U}$, all properties of NMO ellipses examined by Grechka and TsvankinGrechka and Tsvankin (1998) can be derived from the general expressions for ${\bf U}$ given here. Analysis of NMO-velocity surfaces also helps reveal new properties of NMO velocity, which are hidden at a less general level.

Application of these general concepts leads to Dix-type formulae for effective normal-moveout velocity that involve averaging of specific cross-sections of the NMO-velocity surface along the zero-offset ray in heterogeneous anisotropic media. To implement this Dix-type formalism, we derived analytic expressions for computing the NMO-velocity surface and its cross-sections, as well as reconstructing the surface from a single cross-section. It should be emphasized that our averaging procedure does not require any of the slowness components of the zero-offset ray to be preserved between the reflector and the surface. (It is assumed, however, that surfaces of constant slowness, such as boundaries between layers, are locally plane at each point of the ray trajectory.)

Although there exists an alternative way of building the NMO ellipses in heterogeneous media using dynamic ray-tracing equationsGrechka et al. (1999b), the methodology developed here is much more suitable for obtaining closed-form analytic solutions for NMO velocity; it also lends itself to geometrical interpretation. For instance, the Dix-type averaging becomes a purely analytic procedure for the important model composed of homogeneous anisotropic layers separated by plane arbitrary dipping interfaces.

V.2 Implications for anisotropic parameter estimation

The main results of our Dix-type formulation which have important implications in the estimation of anisotropic parameters from reflection data can be summarized as follows:

1. The Dix-type equations operate with the intersections ${\bf W}^{\cal P}$ of the NMO surfaces ${\bf U}$ with generally dipping planes ${\cal P}$ determined by either dipping interfaces along the ray (Figure 4) or the derivative of the slowness vector $d{\bf p}/d\tau_{0}$ (Figure 3).

2. The intersections ${\bf W}^{\cal P}={\bf U}\bigcap{\cal P}$ can depend on the anisotropic parameters that are not constrained by the NMO ellipses in the horizontal plane.

Thus, certain types of lateral heterogeneity may actually help in anisotropic inversion by tilting the planes ${\cal P}$ encountered by the zero-offset ray. For example, non-horizontal cross-sections of the $P$-wave NMO surface in VTI media depend on the individual values of the vertical velocity $V_{P0}$ and the anisotropic parameters $\epsilon$ and $\delta$, while the NMO ellipse in the horizontal plane is controlled by just their combinations $V_{{\rm nmo},P}$ and $\eta$ [equations (18) – (23)]. As a result, the vertical velocity, which determines the depth scale of the model, may be obtained from surface reflection $P$-wave data acquired over a certain class of laterally heterogeneous VTI models. An example presented by Le Stunff et al.Stunff et al. (1999) corroborates this conclusion for a two-layer model containing a VTI layer separated by a dipping interface from an isotropic layer.

Grechka, Pech and TsvankinGrechka et al. (2000a, b) used the Dix-type equations presented here to develop algorithms for $P$-wave stacking-velocity tomography in piecewise-homogeneous VTI media. Their results show that the presence of irregular interfaces may also aid anisotropic parameter estimation in depth (e.g., using reflection tomography) by increasing the angle coverage of reflected rays. A complex subsurface structure, on the other hand, may produce trade-offs between the anisotropic velocity field and the shape of the reflector and intermediate interfaces. Understanding of the properties of NMO-velocity surfaces should help in analyzing these trade-offs and searching for practical ways to overcome them.

VI Conclusions

1. The pure-mode NMO velocity $V_{\rm nmo}$, treated as a function of the direction ${\cal L}$ in 3-D space, forms a quadratic surface that usually is an ellipsoid, an elliptical cylinder, or a one sheeted hyperboloid. If the medium near the common midpoint is homogeneous, the NMO-velocity surface there always has the shape of a cylinder.

2. The NMO ellipse examined by Grechka and TsvankinGrechka and Tsvankin (1998) is the intersection of the NMO-velocity surface with the horizontal plane.

3. The effective NMO ellipse at the surface can be obtained by Dix-type averaging of specifically oriented cross-sections of the NMO-velocity surfaces along the zero-offset ray. This formalism can be applied to any ($P$ or $S$) pure-mode reflection event.

4. The NMO-velocity surface in each anisotropic layer or block encountered by the ray usually depends on more anisotropic parameters than does the intersection of this surface with a horizontal plane (i.e., the NMO ellipse). If the subsurface contains dipping interfaces above the reflector or other types of lateral heterogeneity, these additional parameters contribute to reflection traveltimes measured at the surface and, in some cases, can be estimated using the Dix-type formulae given here. This conclusion is valid for arbitrary anisotropy in each layer, although the parameter-estimation procedure becomes more complicated for lower symmetries.

VII Acknowledgments

We are grateful to members of the A(nisotropy)-Team of the Center for Wave Phenomena (CWP), Colorado School of Mines, for helpful discussions. Careful reviews by Joe Dellinger (BP Amoco), Ken Larner (CSM) and Ray Brown (Oklahoma Geological Survey, associate editor of Geophysics) helped to improve the manuscript. The support for this work was provided by the members of the Consortium Project on Seismic Inverse Methods for Complex Structures at CWP and by the United States Department of Energy (Award #DE-FG03-98ER14908). I. Tsvankin was also supported by the Shell Faculty Career Initiation Grant.

Appendix A Derivation of the NMO-velocity surface

Here we generalize the work of Grechka and TsvankinGrechka and Tsvankin (1998) on the NMO ellipse by developing a small-offset approximation for the squared reflection traveltime $t^{2}$ recorded on a straight common-midpoint (CMP) line parallel to an arbitrary unit vector ${\cal L}$. The derivation is based on expanding the traveltime in a Taylor series in half-offset $h$ near the CMP location ($h=0$). The traveltime field is assumed to be smooth enough for all needed derivatives to exist at zero offset.

If the coordinate of the common midpoint $O$ is denoted by ${\bf y}$ (Figure 6), the coordinates of the source $S$ and receiver $R$ are ${\bf y-x}$ and ${\bf y+x}$, where

[TABLE]

The pure-mode two-way reflection traveltime $t\,$ depends on the positions of the source and receiver and on the coordinate ${\bf r}={\bf r}({\bf y},{\bf x})$ of the reflection point. Summing up the one-way traveltimes $\tau$ corresponding to the downgoing and upgoing rays yields

[TABLE]

Since reflection point dispersal (i.e., the deviation of ${\bf r}$ from ${\bf r}_{0}$ in Figure 6 for $|{\bf x}|\neq 0$) does not change NMO velocityHubral and Krey (1980); Goldin (1986), we assume that nonzero-offset rays are reflected at point ${\bf r}_{0}$,

[TABLE]

Hence, for a given reflection event and a fixed CMP location (constant ${\bf y}$ and ${\bf r}_{0}$), $t$ is a function of the one-way traveltime from the zero-offset reflection point,

[TABLE]

Taking into account equation (31), the first derivative of the traveltime with respect to the half-offset $h$ is given by

[TABLE]

Using equation (33), we find the derivative (34) at zero offset as

[TABLE]

Differentiating equation (34) again yields

[TABLE]

Therefore, at zero offset

[TABLE]

To obtain the equation for NMO velocity along the CMP line ${\cal L}$, we expand the traveltime $t(h)$ in a Taylor series,

[TABLE]

where $t_{0}=2\tau_{0}$ is the two-way zero-offset traveltime. Substituting equations (35) and (37) into equation (38) leads to

[TABLE]

Squaring equation (39) and keeping quadratic and lower-order terms with respect to $h$ yields

[TABLE]

Introducing the source-receiver offset

[TABLE]

we rewrite equation (40) in its final form

[TABLE]

Here the superscript ${\bf T}$ denotes transposition; the $3\times 3$ symmetric matrix ${\bf U}$ is defined as

[TABLE]

and

[TABLE]

are the components of the slowness vector ${\bf p}=[p_{1},p_{2},p_{3}]$.

Comparing equation (42) with the conventional definition of the normal-moveout velocity $V_{\rm nmo}({\cal L})$ on the CMP line ${\cal L}$,

[TABLE]

we conclude that

[TABLE]

Appendix B Constructing the NMO-velocity surface from the NMO ellipse

Let us assume that the NMO ellipse (matrix ${\bf W}$) in the $[x_{1},x_{2}]$-plane has been reconstructed from three or more moveout-velocity measurements in different azimuthal directions. To build the NMO-velocity surface ${\bf U}$ from ${\bf W}$ using equation (5), we need to obtain the matrix elements $U_{k3}$. This can be done by using the Christoffel equation,

[TABLE]

where ${\bf p}={\bf p}({\bf x})$ is the slowness vector of rays generated at the zero-offset reflection point, and $F$ explicitly depends on the spatial coordinates ${\bf x}$ because in heterogeneous media the elastic stiffness coefficients $c_{ij}$ vary in space.

Differentiating equation (47) with respect to $x_{1}$ and $x_{2}$ yields

[TABLE]

where $F_{p_{k}}\equiv\partial F/\partial p_{k}$ and $F_{x_{j}}\equiv\partial F/\partial x_{j}$. The partial derivatives of the horizontal slowness components and the zero-offset traveltime $\tau_{0}$ define the NMO ellipse ${\bf W}$ [see equations (2) and (4)]:

[TABLE]

Substituting equation (49) into equation (48) and solving for $\partial p_{3}/\partial x_{j}$, we find

[TABLE]

Differentiating the Christoffel equation (47) with respect to $x_{3}$ and taking into account that $\partial p_{j}/\partial x_{3}=\partial p_{3}/\partial x_{j}$ because of the symmetry of ${\bf U}$ [equations (2) and (5)] leads to the following expression for $U_{33}$:

[TABLE]

At a fixed spatial location ${\bf x}$ the Christoffel equation can be treated as a relationship between the vertical slowness component $p_{3}\equiv q$ and the horizontal slownesses $p_{1}$ and $p_{2}$. Implicit differentiation of the Christoffel equation then givesGrechka et al. (1999b)

[TABLE]

where $q_{,j}\equiv\partial q/\partial p_{j}\equiv\partial p_{3}/\partial p_{j}$ and $F_{q}\equiv\partial F/\partial q\equiv\partial F/\partial p_{3}$. Using equation (52), we rewrite equations (50) and (51) in the form

[TABLE]

and

[TABLE]

Substituting equation (53) into (54), we obtain the final expression for $U_{33}$ and the NMO-velocity surface as a whole:

[TABLE]

where bullets in the low left-hand corner of the matrix denote the elements $U_{21}=U_{12}$, $U_{31}=U_{13}$ and $U_{32}=U_{23}$.

Appendix C Dix-type averaging in heterogeneous anisotropic media

Here we give a detailed description of the Dix-type procedure for building the NMO-velocity surfaces in heterogeneous anisotropic media. Suppose we would like to use a known surface ${\bf U}(\tau_{0})$ at the zero-offset traveltime $\tau_{0}$ to construct the surface ${\bf U}(\tau_{0}+\Delta\tau_{0})$ at the time $\tau_{0}+\Delta\tau_{0}$, where $\Delta\tau_{0}$ is an infinitesimal interval. We assume that the projection of the slowness vector ${\bf p}(\tau_{0})$ onto the plane ${\cal P}(\tau_{0})\,\bot\,d{\bf p}/d\tau_{0}$ is locally preserved over the ray segment corresponding to $\Delta\tau_{0}$ (Figure 3). If we find the intersection of the NMO-velocity surface ${\bf U}(\tau_{0})$ with the plane ${\cal P}(\tau_{0})$,

[TABLE]

the intersection at the time $\tau_{0}+\Delta\tau_{0}$ is given by the Dix-type equation (27):

[TABLE]

Here

[TABLE]

is the intersection of the *local * NMO-velocity surface {\bf U}\bigl{|}_{[\tau_{0},\,\tau_{0}+\Delta\tau_{0}]} with the plane ${\cal P}(\tau_{0})$. The local surface {\bf U}\bigl{|}_{[\tau_{0},\,\tau_{0}+\Delta\tau_{0}]} can be computed using equations (8) and (15).

The plane {\cal P}(\tau_{0}+\Delta\tau_{0})\,\bot\,d{\bf p}/d\tau_{0}\bigl{|}_{\tau_{0}+\Delta\tau_{0}} at the traveltime $\tau_{0}+\Delta\tau_{0}$ generally differs from the plane ${\cal P}(\tau_{0})$, as shown in Figure 3. In order to account for the rotation of the plane ${\cal P}$ along the ray, we reconstruct the NMO surface ${\bf U}(\tau_{0}+\Delta\tau_{0})$ from its cross-section ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0}+\Delta\tau_{0})$ using equations (65) and (68) (see Appendix D) and find the intersection of ${\bf U}(\tau_{0}+\Delta\tau_{0})$ with the plane ${\cal P}(\tau_{0}+\Delta\tau_{0})$:

[TABLE]

Since the slowness components are locally preserved in the plane ${\cal P}(\tau_{0}+\Delta\tau_{0})$, we can continue the cross-section ${\bf W}^{{\cal P}(\tau_{0}+\Delta\tau_{0})}(\tau_{0}+\Delta\tau_{0})$ along the next ray segment using equation (56) applied to the time interval $[\tau_{0}+\Delta\tau_{0},\,\tau_{0}+2\Delta\tau_{0}]$. Thus, NMO-velocity surfaces in heterogeneous anisotropic media can be computed simultaneously with integrating ray equations (28).

In summary, continuation of the NMO-velocity surface along the zero-offset ray involves the following steps:

Step 1.

Construct the intersection ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0})$ [equation (63) below] of the given NMO-velocity surface ${\bf U}(\tau_{0})$ with the plane ${\cal P}(\tau_{0})$ orthogonal to the vector $d{\bf p}/d\tau_{0}$, which is specified by the second equation (28).

Step 2.

Continue the cross-section of the NMO-velocity surface over the time interval $\Delta\tau_{0}$; that is, compute ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0}+\Delta\tau_{0})$ from ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0})$ using equation (56).

Step 3.

Reconstruct the surface ${\bf U}(\tau_{0}+\Delta\tau_{0})$ [equations (65) and (68) below] from its cross-section ${\bf W}^{{\cal P}(\tau_{0})}(\tau_{0}+\Delta\tau_{0})$.

Step 4.

Repeat Step 1 for the surface ${\bf U}(\tau_{0}+\Delta\tau_{0})$.

Appendix D Operations with cross-sections of NMO-velocity surfaces

In this appendix, we show how to compute the intersection ${\bf W}^{\cal P}$ of the NMO-velocity surface ${\bf U}$ with an arbitrary plane ${\cal P}$ and reconstruct the matrix ${\bf U}$ from its given cross-section ${\bf W}^{\cal P}$. Let us denote by ${\bf z}$ the unit vector in the direction $d{\bf p}/d\tau_{0}$ [see equation (28)] normal to the plane ${\cal P}$. The vector ${\bf z}$ can be specified by two spherical angles $\phi_{1}$ and $\phi_{2}$:

[TABLE]

It is straightforward to verify that the unit vectors

[TABLE]

and

[TABLE]

are both orthogonal to ${\bf z}$ and, therefore, lie in the plane ${\cal P}\perp{\bf z}$. Thus, any vector ${\bf b}$ in ${\cal P}$ is given by

[TABLE]

where $\alpha$ is the azimuth (within ${\cal P}$) with respect to ${\bf b}^{(1)}$.

The NMO velocity [equation (1)] within the plane ${\cal P}$,

[TABLE]

can be viewed as the intersection of the NMO surface ${\bf U}$ with the plane ${\cal P}$. Substituting equations (58) – (60) into (61) yields

[TABLE]

where

[TABLE]

and

[TABLE]

Equations (57) – (59), (63), and (64) define the matrix ${\bf W}^{\cal P}$ that describes the intersection (i.e., the NMO ellipse) of the NMO surface ${\bf U}$ with the plane ${\cal P}$ with the unit normal ${\bf z}$.

Next, we show how to reconstruct the whole NMO surface ${\bf U}$ from its cross-section ${\bf W}^{\cal P}$. This procedure is based on equation (55), which can be written in the form

[TABLE]

where ${\bf W}$ is the NMO ellipse in the horizontal plane $[x_{1},x_{2}]$, and the quantities $A_{i}$ are given by

[TABLE]

[TABLE]

The derivatives of $F$ [equation (47)], which also determine the derivatives $q_{,i}$ and $q_{,ij}$ of the vertical slowness component with respect to the horizontal slownesses [equation (52) and the second equation (15)], are evaluated at point ${\bf x}$ of the zero-offset ray specified by the one-way traveltime $\tau_{0}$.

It is evident from equation (65) that in order to compute ${\bf U}$ we need to know the matrix ${\bf W}$ because all other quantities are obtained by differentiating the Christoffel equation (47). Substituting equation (65) into (63) leads to three linear equations relating the matrices ${\bf W}$ and ${\bf W}^{\cal P}$:

[TABLE]

Here

[TABLE]

and

[TABLE]

To emphasize the fact that equations (66) represent a system of linear equations for the unknown elements $W_{i^{\prime}j^{\prime}}$, we replace the pairs of indexes $\{ij\}$ and $\{i^{\prime}j^{\prime}\}$ by a single index using the following convention: $\{11\}\rightarrow 1$, $\{12\}\rightarrow 2$, and $\{22\}\rightarrow 3$. Then, equations (66) can be rewritten in a more conventional form,

[TABLE]

where $E_{k^{\prime}k}$ are the elements of the $3\times 3$ matrix

[TABLE]

The linear system (67) can be solved for the matrix ${\bf W}$ using standard techniques:

[TABLE]

Thus, equations (65) and (68), supplemented with the Christoffel equation (47), make it possible to reconstruct the NMO-velocity surface ${\bf U}$ from its intersection ${\bf W}^{\cal P}$ with the plane ${\cal P}$.