Efficient depth extrapolation of waves in elastic isotropic media

The Method

We start with the wave equation governing the displacement in an arbitrary heterogenous isotropic elastic medium in the Navier form (see (Segall, 2010)):

$$\rho \ddot{u^i}\;=\;\mu \Delta u^i+ \frac{\mu}{1-2\nu} \frac{ \partial}{\partial x_i} \frac{ \partial u_k}{\partial x_k},\;i=1,2,3,$$ (1)

where $u^i$ denote the components of a displacement field, $\mu$ is the shear modulus, $\nu$ is Poisson's ratio for the medium, and $\rho$ is the density. In this paper we consider a heterogenous elastic medium under the assumption of local homogeneity - otherwise the elastic moduli would not be factored outside of the differentiation operators in equation 1. However, our method can be extended to cover the case when the local homogeneity assumption is dropped. ``Freezing'' the coefficients of equation 1 and applying the Fourier transform in time and horizontal variables $x_1=x,x_2=y$ , and substituting

$$\frac{\mu}{1-2\nu}=\lambda+\mu,$$ (2)

where $\lambda$ is the Lamé coefficient (see (Mavko et al., 2009),(Segall, 2010)), we get

$$\rho \omega^2 u^1+\mu\left[ (-k_x^2-k_y^2)u^1+\frac{\partial^2 u^... ...+\mu)\left[ -k_x^2 u^1-k_x k_y u^2+ik_x\frac{\partial u^3}{\partial z}\right] = 0,$$

$$\rho \omega^2 u^2+\mu\left[ (-k_x^2-k_y^2)u^2+\frac{\partial^2 u^... ...+\mu)\left[ -k_x k_y u^1-k_y^2 u^2+ik_y\frac{\partial u^3}{\partial z}\right] = 0,$$

$$\rho \omega^2 u^3+\mu\left[ (-k_x^2-k_y^2)u^3+\frac{\partial^2 u^... ... \frac{\partial u^2}{\partial z} + \frac{\partial^2 u^3}{\partial z^2}\right] = 0,$$ (3)

where $k_x,k_y$ are horizontal wave numbers and $\omega$ is the frequency. The left-hand side of system 3 is the result of an ordinary differential operator applied to a vector-function $\mathbf{u}=(u^1,u^2,u^3)$ and parameterized by horizontal wave numbers. In the present form equations 3 cannot be used for computationally efficient explicit depth extrapolation in a heterogeneous medium; however, these equations can be used for modeling displacements by solving a series of boundary-value problems (see (Maharramov, 2012)). In (Maharramov, 2012) it was suggested that equations 3 might be factorized in such a way as to allow solving them by alternating one-way extrapolation in opposite directions. More specifically, we seek a factorization of operator equation 3 of the form:

$$\left(E(\lambda,\mu) \frac{\partial}{\partial z} +A(k_x,k_y)+c_{\... ...\mu) \frac{\partial}{\partial z} +B(k_x,k_y)+c_{\omega} I \right) \mathbf{u}=0,$$ (4)

where

$$E(\lambda,\mu)\; = \;\left[ \begin{array}{ccc} \sqrt{\mu} & 0 & 0 \\ 0 & \sqrt{\mu} & 0 \\ 0 & 0 & \sqrt{\lambda+2\mu} \end{array} \right],$$

$$c_{\omega} \; = \; \sqrt{\rho}\omega,$$ (5)

and $A,B$ are $3\times 3$ matrices with components that are complex-valued functions of the horizontal wave numbers, $I$ is the $3\times 3$ identity matrix. Performing the multiplication in equation 4 and using equation 3, we obtain:

$$A(k_x,k_y)B(k_x,k_y)+c_{\omega}[A(k_x,k_y)+B(k_x,k_y)] =P(k_x,k_y),$$

$$A(k_x,k_y) E(\lambda,\mu)+E(\lambda,\mu) B(k_x,k_y)+2c_{\omega} E(\lambda,\mu) =S(k_x,k_y),$$ (6)

where

$$P= \left[ \begin{array}{ccc} -(\lambda+2\mu)k_x^2 - \mu k_y^2 & -(\l... ...+2\mu)k_y^2 - \mu k_x^2 & 0 \\ 0 & 0 & -\mu(k_x^2 + k_y^2) \end{array} \right],$$

$$S= \left[ \begin{array}{ccc} 0 & 0 & i (\lambda+\mu) k_x \\ 0 & 0 & ... ...da+\mu) k_y \\ i(\lambda+\mu) k_x & i(\lambda+\mu) k_y & 0 \end{array} \right].$$ (7)

Combining equations 6 and 7, we get the following equation for the operators $A$ and $B$ :

$$A(k_x,k_y)B(k_x,k_y)+c_{\omega}(A(k_x,k_y)+B(k_x,k_y))=P(k_x,k_y),$$

$$E(\lambda,\mu) B(k_x,k_y)+A(k_x,k_y) E(\lambda,\mu)=\tilde{S}(k_x,k_y),$$ (8)

where

$$\tilde{S}(k_x,k_y)= S(k_x,k_y)-2 c_{\omega} E(\lambda,\mu).$$ (9)

Equations 4, 8 in combination with equations 7 and 9 suggest the following procedure for extrapolating solutions to system 1 in depth:

1. Solve the system of matrix equations 8 for $A,B$ , for each pair of horizontal wave numbers $k_x,k_y$ and two reference values of each elastic parameter $\lambda_{\textrm{min}},\lambda_{\textrm{max}}$ and $\mu_{\textrm{min}},\mu_{\textrm{max}}$ ;
2. Evaluate
   $$\left( E(\lambda,\mu) \frac{\partial}{\partial z}+B(-i\partial_x, -i\partial_y)+ c_{\omega} I \right) \mathbf{u}(x,y,z=0)$$
   from the initial conditions and assign the value to an auxiliary function $\tilde{\mathbf{u}}(x,y,z=0)$ ;
3. Solve

   $$\left( E(\lambda,\mu) \frac{\partial}{\partial z}+A(-i\partial_x, -i\partial_y)+ c_{\omega} I \right) \mathbf{\tilde{u}}(x,y,z)=0$$ (10)

   
   
   by downward continuing in depth, using the formula

   $$\tilde{u}(x,y,z+\Delta z)=\exp{\left[ -\Delta z E^{-1}(A(-i\partial_x,-i \partial_y)+c_{\omega}I)\right]}\mathbf{\tilde{u}}(x,y,z).$$ (11)

   
   

4. Perform each step of the depth extrapolation for four combinations of the reference elastic parameters, then apply the inverse Fourier transform to the four fields and interpolate at each spatial point of the depth slice using true $\lambda(x,y),\mu(x,y)$ as e.g. in the PSPI method (see (Biondi, 2005)).
5. After reaching the desired maximum depth, find the solution $\mathbf{u}$ by upward extrapolation:

   $$\left( E(\lambda,\mu) \frac{\partial}{\partial z} +B(-i\partial_x... ...\partial_y)+ c_{\omega} I \right) \mathbf{u}(x,y,z)= \mathbf{\tilde{u}}(x,y,z).$$ (12)

   
   

6. Repeat the above steps for each frequency component $\mathbf{u}(\omega,x,y,z).$

The above algorithm is stable if the spectrum of matrix

$$A(k_x,k_y)+c_{\omega} I$$ (13)

is not in the interior of the left half-plane, and the spectrum of

$$B(k_x,k_y)+c_{\omega} I$$ (14)

is not in the interior of the right half-plane. While the above algorithm tries to mimic two-way wave propagation, it is effectively just an approximation to the propagation process as it ignores the interaction between the up and down-going wave at intermediate depth steps. A less accurate alternative would be to downward-continue the wave field using equation 11 in a way similar to the one-way depth extrapolation using the scalar square-root equation (see (Claerbout, 1985),(Biondi, 2005)). The latter approach would be unable to image any dips beyond 90$^\circ$ , however, it would reduce the cost of extrapolation by a further factor of 2. Note the cost of solving equation 10 in depth is roughly three times that of solving the scalar square-root equation.

The above analysis may be extended to the case of an arbitrary anisotropic elastic medium. The fact that the components of the pseudo-differential operator matrices $A(-i\partial_x,-i\partial_y), B(-i\partial_x,-i\partial_y)$ are not given in an analytical form, but are only computed numerically, does not limit their applicability.

Factorization of system 3 in the elastostatic case was one of the approaches mentioned by the author in Maharramov (2012). However, the one-way extrapolation technique is mostly useful for elastodynamic problems as the passband of the factorized depth extrapolators (e.g., as in equation 11) narrows down to zero with the temporal frequency passing to the zero static limit.

maximageigenval Figure 1. The phase of a phase-shift operator corresponding to the maximum imaginary part of the eigenvalues of operator 15. Multicomponent ``phase-shift'' is defined by three such scalar phase-shift operators and a $3\times 3$ matrix $Q$ of equation 16.

Note that equation 1 uses elastic parameterization that degenerates into a singularity if the shear modulus is equal to zero. This is not causing any problems with purely acoustic wave extrapolation as the singularity is effectively removed from equations 3 by the substitution in equation 2.

Efficient depth extrapolation of waves in elastic isotropic media