Azimuthally anisotropic 3D velocity continuation

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 (Fomel, 2003b). 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 (1998) 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,

$$T^2(\mathbf{x},\mathbf{y},\tensor{W})=4\left(\tau ^2+\left( ... ...thrm{T}}\tensor{W}\left( \mathbf{x}-\mathbf{y}\right) \right),$$ (1)

where $\tau$ is the one-way vertical traveltime after migration, $\mathbf{x}$ is the $(x_1,x_2)$ surface position of the zero-offset receiver in survey coordinates, $\mathbf{y}$ is the surface position of the point source image, and superscript ${}^{\mathrm{T}}$ denotes transpose. The three independent elements of the symmetric slowness matrix,

$$\tensor{W}=\left( \begin{array}{cc} W_{11} & W_{12} \\ W_{12} & W_{22} \end{array}\right),$$ (2)

have units of slowness-squared, and the eigenvalues and eigenvectors of $\tensor{W}$ determine the symmetry axes of the effective anisotropic medium (Grechka and Tsvankin, 1998). In most common geologic situations, the eigenvalues of $\tensor{W}$ are positive (Tsvankin, 2005), and equation 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. $W_{11}$ and $W_{22}$ are the squared moveout slownesses along their respective survey coordinates, $x_1$ and $x_2$. The third parameter, $W_{12}$, arises from observing the ellipse in the $x_1$-$x_2$ survey coordinates, which are generally rotated relative to its major and minor axes.

The three-parameter moveout model of equation 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 $\tensor{W}$ into more intuitive measurements. If the ellipse happens to be aligned with the survey coordinates, $W_{12}=0$. Finding the rotation angle which properly diagonalizes $\tensor{W}$ therefore allows one to predict the orientation of the symmetry axes. This amounts to an eigenvalue problem, where the fast and slow velocities can be found as the eigenvalues and eigenvectors of $\tensor{W}$. The eigenvalues, $W_{fast}$ and $W_{slow}$, of the slowness matrix can be found following Grechka and Tsvankin (1998),

$$W_{slow,fast}=\frac{1}{2}\left[ W_{11}+W_{22}\pm \sqrt{(W_{11}-W_{22})^2+4W_{12}^2}\right].$$ (3)

Since the eigenvalues have units of slowness squared, the smaller eigenvalue is $W_{fast}=1/v_{fast}^2$. One can solve for the angle $\beta$ between the acquisition coordinates and the symmetry axes by using the formula found by Grechka and Tsvankin (1998),

$$\beta = \tan ^{-1}\left[ \frac{W_{22}-W_{11}+\sqrt{\left(W_{22}-W_{11}\right)^2+4W_{12}^2}}{2W_{12}}\right].$$ (4)

The eigenvalues can then be used together with $\beta$ to solve for the zero-offset migration slowness $S$ as a function of source-receiver azimuth $\theta$:

$$S^2(\theta )=W_{slow}\cos^2(\theta -\beta )+W_{fast}\sin^2(\theta -\beta).$$ (5)

Equations 3-5 allow one to convert the mathematically convenient parameters of $\tensor{W}$ to more intuitive parameters, such as the fastest and slowest propagation velocities ($V_{fast}$,$V_{slow}$), the azimuth of the slowest velocity ($\beta$), and the percent anisotropy ( $\sigma = 100\times \left( 1-V_{slow}/V_{fast}\right)$). Alternatively, $\tensor{W}$ can be converted into other common geophysical parametrizations. For example, Grechka and Tsvankin (1998) show that once the effective parameters $\tensor{W}$ have been converted to slowness as a function of azimuth by equation 5, they can be expressed in terms of horizontal transverse isotropy parameters as,

$$S^2(\theta )=\frac{1}{V_{P0}^2}\frac{1+2\delta ^{(v)}\sin^2(\theta)}{1+2\delta ^{(v)}},$$ (6)

where $\delta ^{(v)}$ is the Thomsen-style parameter (Thomsen, 1986), introduced by Tsvankin (1997), and $V_{P0}$ is the vertical P-wave velocity.

Conventionally, one assumes that equation 1 characterizes a particular event defined in image coordinates ( $\mathbf{x},\tau$), but one can also describe how that event would transform given a change in the image parameters $\tensor{W}$. Regardless of the velocity model, the traveltime $T$ must remain unchanged between different images. From this observation, we arrive at the following set of conditions:

$$\nabla _{\mathbf{x}}T^2 =\left( \begin{array}{c} \frac{\par... ...t)=8\tau \nabla _x\tau +8\tensor{W}(\mathbf{x}-\mathbf{y})=0,$$ (7)

and,

$$\nabla _{\tensor{W}}T^2 = \left( \begin{array}{cc} \frac{\... ...(\mathbf{x}-\mathbf{y})(\mathbf{x}-\mathbf{y})^{\mathrm{T}}=0.$$ (8)

Combining and reducing these conditions yields a system of equations that are defined only in the image parameter coordinates,

$$2 \frac{\partial \tau }{\partial W_{11}}+\frac{\tau \left(W_... ...artial x_2}\right)^2}{\left(W_{12}^2-W_{11}W_{22}\right)^2}=0,$$ (9)

$$2 \frac{\partial \tau }{\partial W_{22}}+\frac{\tau\left(W_{... ...artial x_2}\right)^2}{\left(W_{12}^2-W_{11}W_{22}\right)^2}=0,$$ (10)

and,

$$2 \frac{\partial \tau }{\partial W_{12}}-\frac{2\tau \left(W... ...\partial x_2}\right)}{\left(W_{12}^2-W_{11}W_{22}\right)^2}=0.$$ (11)

The system of kinematic equations describing azimuthally anisotropic velocity continuation is then found by combining equations 9-11. In a vector notation, this becomes

$$\nabla _{\tensor{W}}\tau +\frac{\tau }{2}\tensor{W}^{-1} \na... ...nabla _{\mathbf{x}}\tau \right)^{\mathrm{T}}\tensor{W}^{-1}=0,$$ (12)

where $\nabla _{\mathbf{x}}$ and $\nabla _{\tensor{W}}$ are in the form given by equations 7 and 8.

The method of characteristics (Courant and Hilbert, 1989) provides a link between a kinematic equation (such as 12) and its corresponding wave-type equation. Fomel (2003b) 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,

$$\psi =t-\tau (\mathbf{x},\tensor{W})=0,$$ (13)

and replacing $\tau$ with $\psi$ and $t$, we obtain an alternative form of equation 12,

$$\psi _t\nabla _{\tensor{W}}\psi +\frac{t}{2}\tensor{W}^{-1} ... ...nabla _{\mathbf{x}}\psi \right)^{\mathrm{T}}\tensor{W}^{-1}=0.$$ (14)

Equation 13 guarantees that the wavefronts of time-domain image wavefield $P$ exist only where the arrival time $\tau$ is equal to the recorded time $t$ at a given location. Now take both $\xi _i$ and $\xi _j$ to represent each of $t$, $W_{11}$, $W_{12}$, $W_{22}$, $x_1$, and $x_2$. According to the method of characteristics, if $\Lambda _{ij}$ is the coefficient in front of $\frac{\partial \psi }{\partial \xi _i}\frac{\partial \psi }{\partial \xi _j}$ from kinematic equation 14, then the corresponding wave equation will have the same coefficients $\Lambda _{ij}$ in front of each $\frac{\partial ^2P }{\partial \xi _i\partial \xi _j}$ derivative. The time-derivative $\psi _t$ is equal to $1$ given equation 13, and is included in the first term of equation 14 to facilitate the use of the method of characteristics. Then, by introducing $\tensor{P}_{xx}$ as the spatial Hessian matrix of the wavefield,

$$\tensor{P}_{xx}=\left( \begin{array}{cc} \frac{\partial ^2P... ...1} & \frac{\partial ^2P }{\partial x_2^2} \end{array}\right),$$ (15)

we arrive at the azimuthally anisotropic post-stack velocity continuation wave equation,

$$\nabla _{\tensor{W}}P_t=-\frac{t}{2}\tensor{W}^{-1}\tensor{P}_{xx}\tensor{W}^{-1}.$$ (16)

In the isotropic case, $\tensor{W}$ is diagonal and $W_{11}=W_{22}$. Equation 16 then reduces to the isotropic velocity continuation equation first derived by Claerbout (1986).

Azimuthally anisotropic 3D velocity continuation