Recursive integral time extrapolation of elastic waves using low-rank symbol approximation

Homogeneous Media

We first assume homogeneous model properties and investigate the form of time stepping solutions in general elastic media. For elastic anisotropic wave equations, the matrix operator $\mathbf{A}$ corresponds to the spatial Fourier transform of the matrix $-\Gamma /\rho = -\mathbf{D}\mathbf{C}\mathbf{D}^\intercal/\rho$ , where $\Gamma$ is the Christoffel matrix, $\rho$ is density, $\mathbf{C}$ is the elastic stiffness tensor expressed in Voigt notation as a $6\times6$ matrix, and $\mathbf{D}$ is the derivative matrix operator give by

$$\mathbf{D}= \begin{bmatrix}\partial_x & 0 & 0 & 0 & \partial_z & ... ...partial_x \\ 0 & 0 & \partial_z & \partial_y & \partial_x & 0 \end{bmatrix} \;.$$ (11)

In the Fourier (wavenumber) domain, the homogeneous $\mathbf{A}$ takes the form $\tilde{\Gamma }/\rho$ after $i^2$ cancels the negative sign, i.e., $-\Gamma \xrightarrow{\mathscr{F}} \tilde{\Gamma }$ . For example, the Christoffel matrix $\tilde{\Gamma }$ in the case of orthorhombic anisotropy takes the form:

$$\tilde{\Gamma } = \begin{bmatrix}C_{11}k_x^2+C_{66}k_y^2+C_{55}k_... ...& (C_{23}+C_{44})k_yk_z & C_{55}k_x^2+C_{44}k_y^2+C_{33}k_z^2 \end{bmatrix} \;.$$ (12)

The square root matrix $\Phi$ is analogous to the phase function in the acoustic case, and corresponds to the angular frequency $\omega$ according to the dispersion relation. Since the matrix $\tilde{\Gamma }$ is symmetric positive definite (SPD), it can be diagonalized with its eigenvalues corresponding to the square of phase velocity of separate wave modes and its orthogonal eigenvectors corresponding to the polarization directions:

$$\mathbf{A}= \frac{\tilde{\Gamma }}{\rho} = \mathbf{Q} \mathbf{V} \mathbf{Q}^\intercal$$ (13)

$$= \begin{bmatrix}\mathbf{a}_p & \mathbf{a}_{s1} & \mathbf{a}_{s2}\e... ...cal \\ \mathbf{a}_{s1}^\intercal \\ \mathbf{a}_{s2}^\intercal \end{bmatrix} \;.$$

Since $\mathbf{Q}$ is orthogonal, $\mathbf{Q}^\intercal$ projects the input vector to its column space. The square root of $\tilde{\Gamma }$ is found by taking the square root of the eigenvalues in the diagonal matrix. Analogously, the wave extrapolation operator, $e^{i\Phi \Delta t}$ , can be computed as:

$$e^{i\Phi \Delta t} =$$ (14)

$$\begin{bmatrix}\mathbf{a}_p & \mathbf{a}_{s1} & \mathbf{a}_{s2}\e... ...cal \\ \mathbf{a}_{s1}^\intercal \\ \mathbf{a}_{s2}^\intercal \end{bmatrix} \;.$$

An analogous idea of using matrix exponentials for wave propagation was studied previously by Kosloff and Kessler (1987) in application to the one-way wave equation. Physically, the operator in equation 14 first decomposes the input vector wavefield into three wave modes, phase shifts them using the corresponding phase velocities, and then aligns them in the polarization directions of the decomposed wave modes. The operator defined in equation 14 can be expressed as a summation of rank-one matrices

$$e^{i\Phi \Delta t} = \sum\limits_{i=p,s1,s2} e^{i v_i k \Delta t} \mathbf{a}_i \mathbf{a}_i^\intercal \; .$$ (15)

Note that $\mathbf{a}_i \mathbf{a}_i^\intercal$ term in equation 15 is the wave mode decomposition operator (Zhang and McMechan, 2010; Dellinger, 1991).

In 3D, we can transform the input vector in the wavenumber domain $\hat{\mathbf{u}}(k) = \begin{bmatrix}\hat{u}_x(k) & \hat{u}_y(k) & \hat{u}_z(k) \end{bmatrix}^\intercal$ . To apply the Fourier integral operator, we omit the pair of forward and backward Fourier transforms and formally write

$$e^{i\Phi \Delta t} \hat{\mathbf{u}}(k) = \sum\limits_{i=p,s1,s2} e^{i v_i k \Delta t} \mathbf{a}_i \mathbf{a}_i^\intercal \hat{\mathbf{u}}$$ (16)

$$= \begin{bmatrix}\mathbf{a}_p & \mathbf{a}_{s1} & \mathbf{a}_{s2}\e... ...f{a}_{s1}^\intercal \\ \mathbf{a}_{s2}^\intercal \end{bmatrix} \hat{\mathbf{u}}$$

$$= \begin{bmatrix} s_{xx} & s_{xy} & s_{xz} \\ s_{yx} & s_{yy} & s_{yz} \\ s_{zx} & s_{zy} & s_{zz} \\ \end{bmatrix} \hat{\mathbf{u}} \;,$$

where

$$s_{xx} = e^{i v_p k \Delta t} a_{p_x} a_{p_x} + e^{i v_{s1} k \Delta t} a_{s1_x} a_{s1_x} + e^{i v_{s2} k \Delta t} a_{s2_x} a_{s2_x} \;,$$ (17)

$$s_{xy} = e^{i v_p k \Delta t} a_{p_x} a_{p_y} + e^{i v_{s1} k \Delta t} a_{s1_x} a_{s1_y} + e^{i v_{s2} k \Delta t} a_{s2_x} a_{s2_y} \;,$$

$$s_{xz} = e^{i v_p k \Delta t} a_{p_x} a_{p_z} + e^{i v_{s1} k \Delta t} a_{s1_x} a_{s1_z} + e^{i v_{s2} k \Delta t} a_{s2_x} a_{s2_z} \;,$$

$$s_{yy} = e^{i v_p k \Delta t} a_{p_y} a_{p_y} + e^{i v_{s1} k \Delta t} a_{s1_y} a_{s1_y} + e^{i v_{s2} k \Delta t} a_{s2_y} a_{s2_y} \;,$$

$$s_{yz} = e^{i v_p k \Delta t} a_{p_y} a_{p_z} + e^{i v_{s1} k \Delta t} a_{s1_y} a_{s1_z} + e^{i v_{s2} k \Delta t} a_{s2_y} a_{s2_z} \;,$$

$$s_{zz} = e^{i v_p k \Delta t} a_{p_z} a_{p_z} + e^{i v_{s1} k \Delta t} a_{s1_z} a_{s1_z} + e^{i v_{s2} k \Delta t} a_{s2_z} a_{s2_z} \;,$$

$$s_{yx} = s_{xy} \;,\; s_{zx} = s_{xz} \;,\; s_{zy} = s_{yz} \;.$$

Using the form of equation 9, we can also define the backward time extrapolator

$$\hat{\mathbf{u}}(\mathbf{x},t-\Delta t) = e^{-i\Phi \Delta t} \hat{\mathbf{u}}(\mathbf{x},t) \;.$$ (18)

If we require the data to be real-valued at every time step, and sum the forward and backward extrapolators, we arrive at the two-step formulation:

$$\hat{\mathbf{u}} (\mathbf{x},t+\Delta t) = 2\,\cos(\Phi \Delta t) \hat{\mathbf{u}}(\mathbf{x},t) - \hat{\mathbf{u}}(\mathbf{x},t-\Delta t) \;.$$ (19)

According to equation 15, the cosine term has the following interpretation:

$$\cos(\Phi \Delta t) = e^{i\Phi \Delta t} + e^{-i\Phi \Delta t}$$ (20)

$$= \sum\limits_{i=p,s1,s2} (e^{i v_i k \Delta t} + e^{-i v_i k \Delta t}) \mathbf{a}_i \mathbf{a}_i^\intercal$$

$$= \sum\limits_{i=p,s1,s2} 2\,\cos(v_i k \Delta t) \mathbf{a}_i \mathbf{a}_i^\intercal \;,$$

and can be calculated as

$$\cos(\Phi \Delta t) \hat{\mathbf{u}}(k) = \begin{bmatrix} c_{xx} & c_{xy} & c_{xz} \\ c_{yx} & c_{yy} & c_{yz} \\ c_{zx} & c_{zy} & c_{zz} \\ \end{bmatrix} \hat{\mathbf{u}} \;.$$ (21)

where

$$c_{xx} = \cos(v_p k \Delta t) a_{p_x} a_{p_x} + \cos(v_{s1} k \Delta t) a_{s1_x} a_{s1_x} + \cos(v_{s2} k \Delta t) a_{s2_x} a_{s2_x} \;,$$ (22)

$$c_{xy} = \cos(v_p k \Delta t) a_{p_x} a_{p_y} + \cos(v_{s1} k \Delta t) a_{s1_x} a_{s1_y} + \cos(v_{s2} k \Delta t) a_{s2_x} a_{s2_y} \;,$$

$$c_{xz} = \cos(v_p k \Delta t) a_{p_x} a_{p_z} + \cos(v_{s1} k \Delta t) a_{s1_x} a_{s1_z} + \cos(v_{s2} k \Delta t) a_{s2_x} a_{s2_z} \;,$$

$$c_{yy} = \cos(v_p k \Delta t) a_{p_y} a_{p_y} + \cos(v_{s1} k \Delta t) a_{s1_y} a_{s1_y} + \cos(v_{s2} k \Delta t) a_{s2_y} a_{s2_y} \;,$$

$$c_{yz} = \cos(v_p k \Delta t) a_{p_y} a_{p_z} + \cos(v_{s1} k \Delta t) a_{s1_y} a_{s1_z} + \cos(v_{s2} k \Delta t) a_{s2_y} a_{s2_z} \;,$$

$$c_{zz} = \cos(v_p k \Delta t) a_{p_z} a_{p_z} + \cos(v_{s1} k \Delta t) a_{s1_z} a_{s1_z} + \cos(v_{s2} k \Delta t) a_{s2_z} a_{s2_z} \;,$$

$$c_{yx} = c_{xy} \;,\; c_{zx} = c_{xz} \;,\; c_{zy} = c_{yz} \;.$$

Note that, because the elements of the $3\times3$ matrix $\cos(\Phi \Delta t)$ in equation 21 are real-valued, the wavefield can stay real-valued as well. The matrix $\cos(\Phi \Delta t)$ is also closely connected with the k-space adjustment of the Christoffel matrix (Liu, 1995; Firouzi et al., 2012; Cheng et al., 2016).

Expanding the wave extrapolation operator in the form of equations 17 and 22 reveals a simple relationship between wave propagation and wave mode decomposition. In the one-step formulation, each individual term contained in each element of $e^{i\Phi \Delta t}$ is essentially a sequence of wave-mode decomposition, phase shift and recomposition. For example, the action of the first term in $s_{xx}$ , $e^{i v_p k \Delta t} a_{p_x} a_{p_x}$ , can be interpreted as projecting the x-component of a vector elastic wavefield onto the P-wave mode, phase shifting it using the P-wave phase velocity and then aligning it with the x-component of the P-wave polarization direction. If individual wave modes are needed to perform imaging, the decoupled operators can be separately applied. The operators concerning a specific wave mode are indicated by the corresponding phase velocity. The first columns in equations 17 and 22 are P-wave propagators, while the second and third columns are $S1$ and $S2$ wave propagators, respectively. In general anisotropic media beyond tilted transverse isotropic (TTI) symmetry, the two S-wave modes do not decouple easily. Therefore, in order to avoid S-wave singularities in wave propagation, the two coupled S-waves should be computed together (Sun et al., 2016a; Cheng et al., 2016).

The proposed framework, however, does not require explicit wave mode decomposition for wave extrapolation. Therefore, it can significantly reduce the computational cost of wave extrapolation, yet still obtain waves free of instability and dispersion artifacts. We also emphasize that the proposed method is capable of handling general anisotropic media, including the case of triclinic anisotropy.

Recursive integral time extrapolation of elastic waves using low-rank symbol approximation