Fast algorithms for elastic-wave-mode separation and vector decomposition using low-rank approximation for anisotropic media

Low-rank approximation solutions

We observe that both elastic wave mode separation and vector decomposition are based on polarization of wave modes, and any wave mode shares the algorithm structure of separation or decomposition. We will use qP-wave as an example.

We first write the equivalent version of the first equation of equation 4 in the space-domain as

$$qP(\mathbf{x})=\int{e^{i\mathbf{k}\mathbf{x}} \left[ia_{px}(\math... ...})+ia_{pz}(\mathbf{k})\tilde{U}_{z}(\mathbf{k})\right] }\,\mathrm{d}\mathbf{k}.$$ (12)

To tackle spatial variations of the polarization in heterogeneous media, we extend the integral operators using

$$qP(\mathbf{x}) = \int{ e^{i\mathbf{k}\mathbf{x}}\left[ia_{px}(\mathbf{x},\mathbf{k... ...\mathbf{x},\mathbf{k})\right] \tilde{U}_{y}(\mathbf{k}) }\,\mathrm{d}\mathbf{k}$$

$$+ \int{ e^{i\mathbf{k}\mathbf{x}}\left[ia_{pz}(\mathbf{x},\mathbf{k})\right] \tilde{U}_{z}(\mathbf{k}) }\,\mathrm{d}\mathbf{k},$$ (13)

where $a_{px}(\mathbf {x},\mathbf {k})$ , $a_{py}(\mathbf{x},\mathbf{k})$ and $a_{pz}(\mathbf {x},\mathbf {k})$ represent the $x$ -, $y$ - and $z$ -components of the normalized polarization vector of qP waves at location $\mathbf{x}$ . Compared with nonstationary filtering (Yan and Sava, 2009b), these pseudospectral-like operations are more accurate but less efficient. The computation complexity of the straight-forward implementation is $O(N_{\mathbf{x}}^2)$ , which is prohibitively expensive when the size of model $N_{\mathbf{x}}$ is large.

Similarly, from equation 9, we can derive the space-wavenumber-domain operators for decomposing qP-waves. For example, the $x$ -component of qP-wave satisfies,

$$U_{x}^{qP}(\mathbf{x}) = \int{e^{i\mathbf{k}\mathbf{x}}\left[a_{px}(\mathbf{x},\mathbf{k})... ...}(\mathbf{x},\mathbf{k})\right]\tilde{U}_{y}(\mathbf{k})}\,\mathrm{d}\mathbf{k}$$

$$+ \int{ e^{i\mathbf{k}\mathbf{x}}\left[a_{px}(\mathbf{x},\mathbf{k}... ...\mathbf{x},\mathbf{k})\right] \tilde{U}_{z}(\mathbf{k})}\,\mathrm{d}\mathbf{k}.$$ (14)

Note that more multiplication operations are needed for vector decomposition.

The discrete implementation of each integral operation in equation 13 or 14 naturally arises as a numerical approximation of a continuous Fourier integral operator (FIO) of the general form. Underlying fast solutions to FIOs is a mathematical insight concerning restriction of the integral kernel to subsets of space and wavenumber domains. Whenever these subsets obey a simple geometric condition, the restricted kernel is approximately low rank (Candes et al., 2007; Ying, 2012). Recently, several two-way wave extrapolation operators have been developed with the help of a low-rank approximation of the space-wavenumber-domain wave-propagator matrix in variable and possibly anisotropic media (Alkhalifah, 2013; Fomel et al., 2010; Song and Alkhalifah, 2013; Song et al., 2011; Fomel et al., 2013).

For both mode separation and vector decomposition, phase terms in the FIOs are relatively simple and can be absorbed into inverse Fourier transforms. Therefore our main task is to respectively construct low-rank decomposition for the amplitude terms in the kernel. Previously, Fomel et al. (2013) applied low-rank approximation to the phase-only terms in wave extrapolation operators. Any amplitude term bracketed in equations 13 and 14 can be approximated by the following separated representation:

$$W_{j}(\mathbf{x},\mathbf{k})\approx \sum_{m=1}^M \sum_{n=1}^N B(\mathbf{x},\mathbf{k}_{m})A_{mn}C(\mathbf{x}_{n},\mathbf{k}),$$ (15)

in which $j=x, y, z$ , $N$ and $M$ represent the rank of this decomposition, $W_{j}(\mathbf{x},\mathbf{k})$ represents the mixed-domain mode separation or vector decomposition operator, $B(\mathbf{x},\mathbf{k}_{m})$ is a mixed-domain matrix with reduced wavenumber dimension, $C(\mathbf{x}_{n},\mathbf{k})$ is a mixed-domain matrix with reduced spatial dimension, and $A_{mn}$ is a ${M}\times{N}$ matrix. The construction of the separated form 15 follows the method of Engquist and Ying (2009). It can be viewed as a matrix decomposition problem, i.e.,

$$\mathbf{W} \approx \mathbf{B}\mathbf{A}\mathbf{C},$$ (16)

where $\mathbf{W}$ is the $N_{\mathbf{x}}\times{N_{\mathbf{x}}}$ matrix with entries $W_{j}(\mathbf{x},\mathbf{k})$ , $\mathbf{B}$ is the submatrix of $\mathbf{W}$ that consists of columns associated with $\{\mathbf{k}_{m}\}$ , $\mathbf{C}$ is the submatrix that consists of rows associated with $\{\mathbf{x}_{n}\}$ , and $\mathbf{A}=\{A_{mn}\}$ . Physically, a separable low-rank approximation amounts to selecting a set of $N$ ( $N\ll{N_{x}}$ ) representative spatial locations and $M$ ( $M\ll{N_{x}}$ ) representative wavenumbers. As explained by Fomel et al. (2013) in detail, we first need to restrict the mixed-domain $\mathbf{W}$ to $n$ randomly selected rows. In practice, $n$ can be scaled as $O(r\log{N_{x}})$ and $r$ represents the numerical rank of $\mathbf{W}$ . Then we apply pivoted QR algorithm (Golub and Loan, 1996) to find the corresponding columns for $B(\mathbf{x},\mathbf{k}_{m})$ . To find the rows for $C(\mathbf{x}_{n},\mathbf{k})$ , we apply the pivoted QR algorithm to $\mathbf{\mathbf{W}^{\ast}}$ . The algorithm does not require, at any step, access to the full-matrix $\mathbf{W}$ , only to its selected rows and columns.

Representation 15 speeds up the computation of the FIOs in equations 13 and 14 since

$$\begin{split} &\int{e^{i\mathbf{k}\mathbf{x}}W_{j}(\mathbf{x}... ...(\mathbf{k})\,\mathrm{d}\mathbf{k}}\right)\right). \end{split}$$ (17)

The evaluation of the last formula is effectively equivalent to applying $N$ inverse fast Fourier transforms (FFTs). Accordingly, with low-rank approximation, the computation complexity reduces to $O(NN_{\mathbf{x}}\log{N_{\mathbf{x}}})$ . In other word, the costs are mainly controled by the model size $N_{x}$ and the rank $N$ , which depends on the complexity of the anisotropic velocity model. For isotropic models with arbitrary heterogeneity, the rank automatically reduces to $1$ because the polarization directions are material-independent. Similarly to the observation by Fomel et al. (2013), there is a natural tradeoff in the selection of N: larger values lead to a more accurate separated representation but require a longer computational time. In the examples of the next section, the ranks are automatically calculated based on the estimate of the approximation accuracy and generally aiming for the relative single-precision accuracy (namely the maximum allowable error in low-rank decomposition) of $10^{-6}.$ In multiple-core implementations, the matrix operations in equation 17 are easy to parallelize.

For some applications such as ERTM in TI media, one should construct the separated representations of the operator matrixes in advance, and then implement mode separation or/and decomposition of the extrapolated elastic wavefields before applying the imaging condition. To further save computational cost, appropriate relaxing of the accuracy requirement for low-rank approximation and applying mode separation only every two or three time steps are both good choices in practice.

Fast algorithms for elastic-wave-mode separation and vector decomposition using low-rank approximation for anisotropic media