Simulating propagation of separated wave modes in general anisotropic media, Part I: qP-wave propagators

CORRECTION OF PROJECTION DEVIATION TO REMOVE qS-WAVES

According to the theory of wave mode separation in anisotropic media, one needs to project the elastic wavefields onto the polarization direction to get the separated wavefields of the given mode (Dellinger and Etgen, 1990). Mathematically, this can be implemented through a dot product of the original vector wavefields and the polarization vector in the wavenumber domain (Dellinger and Etgen, 1990) or applying pseudo-derivative operators to the vector wavefields in the space-domain (Yan and Sava, 2009). However, the pseudo-pure-mode qP-wave equations are derived by a similarity transformation aiming to project the displacement wavefield onto the isotropic reference of qP-wave's polarization vector. A partial mode separation has been automatically achieved during wavefield extrapolation using the pseudo-pure-mode qP-wave equations. For typical anisotropic earth media, thanks to the small departure of qP-wave's polarization direction from its isotropic reference, the resulting pseudo-pure-mode wavefields are dominated by qP-wave energy and contaminated by residual qS-waves. To achieve a complete mode separation, we should further correct the projection deviations resulting from the differences between polarization and its isotropic reference. In other words, we split the conventional one-step wave mode separation for anisotropic media (Yan and Sava, 2009; Dellinger and Etgen, 1990) into two steps, of which the first one is implicitly achieved during extrapolating the pseudo-pure-mode wavefields and the second one is implemented after that using the approach that we will present immediately.

Taking VTI as an example, the deviation angle $\zeta$ between the polarization and propagation directions has a complicated nonlinear relation with anisotropic parameters and the phase angle (see Appendix C). According to its expression for weak anisotropic VTI media (Rommel, 1994; Tsvankin, 2001), it seems that the deviation is mainly affected by the difference between $\epsilon$ and $\delta$ , the magnitude of $\delta$ (when $\epsilon-\delta$ stays the same) and the ratio of vertical velocities of qP- and qS-wave, as well as the phase angle. It is possible to design a filtering algorithm to suppress the residual qS-waves using the deviation angle given under the assumption of weak anisotropy. To completely remove the residual qS-waves and correctly separate the qP-waves for arbitrary anisotropy, we propose an accurate correction approach according to the deviation between polarization and wave vectors.

Considering equations 8, 9, 13, and 15, we first decompose the polarization vector of qP-wave $\mathbf {a}_{p}$ as follows:

$$\mathbf{a}_{p}=\mathbf{E}_{p}\mathbf{k},$$ (29)

where the deviation operator satisfies,

$$\mathbf{E}_{p}= \begin{pmatrix}\frac{a_{px}}{k_{x}} &0 &0 \cr 0 & \frac{a_{py}}{k_{y}} &0 \cr 0 & 0 & \frac{a_{pz}}{k_{z}}\end{pmatrix}.$$ (30)

This matrix can be constructed once the qP-wave polarization directions are determined based on the local medium properties at a grid point. For TI media, there are analytical expressions for the qP-wave polarization vectors (Dellinger, 1991). For other anisotropic media with lower symmetry (such as orthorhombic media), we have to numerically compute the polarization vectors using the Christoffel equation.

Then we correct the pseudo-pure-mode qP-wave fields $\widetilde{\overline{\mathbf{u}}}=(\widetilde{\overline{u}}_{x}, \widetilde{\overline{u}}_{y}, \widetilde{\overline{u}}_{z})^{T}$ using a wavenumber-domain filtering based on the deviation operator:

$$\widetilde{\mathbf{u}}_{p}=\mathbf{E}_{p}\cdot{\widetilde{\overline{\mathbf{u}}}},$$ (31)

and finally extract the scalar qP-wave data using

$$u_{p}=u_{px}+u_{py}+u_{pz}$$ (32)

after 3D inverse Fourier transforms. Here, the magnitude of the deviation operator for a certain wavenumber $k=\sqrt{k^2_{x}+k^2_{y}+k^2_{z}}$ is a constant becuase this operator is computed by using the normalized wave and polarization vectors in equation 30. This ensures that for a certain wavenumber, the separated qP-waves are uniformly scaled. More important, this correction step thoroughly removes the residual qS-wave energy.

In heterogeneous anisotropic media, the polarization directions and thus the deviation operators vary spatially, depending on the local material parameters. To account for spatial variability, we propose an equivalent expression to equations 31 and 32 as a nonstationary filtering in the space domain at each location,

$$u_{p}=E_{px}(\overline{u}_{x})+E_{py}(\overline{u}_{y})+E_{pz}(\overline{u}_{z})$$ (33)

where the pseudo-derivative operators $E_{px}(\cdot)$ , $E_{py}(\cdot)$ , and $E_{pz}(\cdot)$ represent the inverse Fourier transforms of the diagonal elements in the deviation matrix $\mathbf {E}_{p}$ .

Figure 1 displays the wavenumber-domain operators of projection onto isotropic (reference) and anisotropic polarization vectors (namely $\mathbf {k}$ and $\mathbf {a}_{p}$ ) as well as the corresponding deviation operator $\mathbf {E}_{p}$ for a 2D homogeneous VTI medium with $v_{p0}=3000m/s$ , $v_{s0}=1500m/s$ , $\epsilon=0.25$ and $\delta=-0.25$ . Note that $\mathbf {E}_{p}$ is not simply the difference between $\mathbf {k}$ and $\mathbf {a}_{p}$ , and $\mathbf {E}_{p}$ becomes the identity operator in case of an isotropic medium. In the space-domain, projecting onto isotropic polarization directions is equivalent to a divergence operation using partial derivative operators, while projection onto polarization directions of qP-waves use operators that have the character of pseudo-derivative operators, due to anisotropy (see Figure 2). Figure 3 shows that the variation of the anisotropy changes the deviation operators greatly. The weaker the anisotropy, the more compact the deviation operators appear. The observation is basically consistent with the equation of polarization deviation angle for VTI media with weak anisotropy. The exact pseudo-derivative operators are very long series in the discretized space domain. Generally, the far ends of these operators have ignorable values even for strong anisotropy. Therefore, in practice, we could truncate the operators to make the spatial filters short and computationally efficient.

This procedure to separate qP-waves, although accurate, is computationally expensive, especially in 3D heterogeneous media. Like the computational problem in conventional wave mode separation from the anisotropic elastic wavefields (Yan and Sava, 2009,2012), the spatial filtering to separate qP-waves is significantly more expensive than extrapolating the pseudo-pure-mode wavefields. In practice, we find that it is not necessary to apply the filtering at every time step. A larger time interval is allowed to save costs enormously, especially for RTM of multi-shot seismic data. According to our experiments, there is little difference between the two migrated images when the filtering is applied at every one and two time step (if only the filtered wavefields are used in the imaging procedure), although about three-forth of the original computational cost are saved for the filtering in the latter case. Moreover, filtering at every three time step still produces an acceptable migrated image. We may further improve the efficiency of the filtering procedure by using the algorithm that resembles the phase-shift plus interpolation (PSPI) scheme recently used in anisotropic wave mode separation (Yan and Sava, 2011). Alternatively, we may greatly reduce the compuational cost but guarantee the accuracy using the mixed (space-wavenumber) domain filtering algorithm based on low-rank approximation (Cheng and Fomel, 2013).

adxNT,apxNT,apvxNT,adzNT,apzNT,apvzNT
Figure 1. Normalized wavenumber-domain operators of projection onto isotropic (reference) and anisotropic polarization vectors of qP-waves, and wavenumber-domain deviation operators in a 2D homogeneous VTI medium: $\mathbf {k}$ (left), $\mathbf {a}_{p}$ (middle) and $\mathbf {E}_{p}$ (right); Top: x-component, Bottom: z-component.

adxx,apxx,apvxx,adzz,apzz,apvzz
Figure 2. Space-domain operators of projecting onto isotropic (left) and anisotropic (middle) polarization vectors, and the corresponding deviation operators (right): Top: x-component, Bottom: z-component. Note that the operators are tapered before transforming into space-domain and the same gain is applied to all pictures to highlight the differences among these operators.

apvxx1,apvxx2,apvxx3,apvxx4,apvxx5,apvzz1,apvzz2,apvzz3,apvzz4,apvzz5
Figure 3. Comparison of the spatial domain deviation operators in VTI media with varied anisotropy strength: In all cases, $v_{p0}=3000m/s$ , $v_{s0}=1500m/s$ , and $\epsilon$ is fixed as 0.2. From left to right, $\delta$ is set as 0.2, 0.1, 0, -0.1, and -0.2, respectively. Top: x-components; Bottom: z-components. To highlight the differences,the same gain is applied to all pictures.

In kinematics, it seems that we can extract scalar qSV-wave fields from the pseudo-pure-mode qP-wave fields $\overline{\mathbf{u}}$ by filtering according to the projection deviation defined by qSV-wave's polarization and wave vector. However, unlike separation of the qP-wave mode, the large projection deviations for qSV-wave modes would result in significant discontinuities in the wavenumber-domain correction operators and strong tails extending off to infinity in the space domain. Accordingly, this reduces compactness of the spatial filters, which prohibits applying the same truncation as for qP-wave spatial filters to reduce computational cost. Computational tricks such as smoothing may result in distorted and imcomplete separation. That is why we are developing a similar approach to simulate propagation of separated qS-wave modes based on their own pseudo-pure-mode wave equations and the corresponding projection deviation corrections for anisotropic media (Kang and Cheng, 2012).

Simulating propagation of separated wave modes in general anisotropic media, Part I: qP-wave propagators