Efficient anisotropic quasi-P wavefield extrapolation using an isotropic low-rank approximation

Efficient anisotropic quasi-P wavefield extrapolation using an isotropic low-rank approximation Summary The computational cost of quasi-P wave extrapolation depends on the complexity of the medium, and specifically the anisotropy. Our effective-model method splits the anisotropic dispersion relation into an isotropic background and a correction factor to handle this dependency. The correction term depends on the slope (measured using the gradient) of current wavefields and the anisotropy. As a result, the computational cost is independent of the nature of anisotropy, which makes the extrapolation efficient. A dynamic implementation of this approach decomposes the original pseudo-differential operator into a Laplacian, handled using the low-rank approximation of the spectral operator, plus an angular dependent correction factor applied in the space domain to correct for anisotropy. We analyse the role played by the correction factor and propose a new spherical decomposition of the dispersion relation. The proposed method provides accurate wavefields in phase and more balanced amplitudes than a previous spherical decomposition. Also, it is free of SV-wave artefacts. Applications to a simple homogeneous transverse isotropic medium with a vertical symmetry axis (VTI) and a modified Hess VTI model demonstrate the effectiveness of the approach. The Reverse Time Migration applied to a modified BP VTI model reveals that the anisotropic migration using the proposed modelling engine performs better than an isotropic migration. Numerical modelling, Seismic anisotropy, Wave propagation 1 INTRODUCTION Seismic anisotropy has been recognized as an important factor in seismic imaging and inversion. Because of the limitations we face in our computational ability and those we face in estimating anisotropy, the acoustic approximation (Alkhalifah 1998) is still the most efficient and practical way to incorporate anisotropy into regular P-wave processing routines. Compared to solving the elastic wave equations, solving the acoustic anisotropic wave equations has two main advantages, which are less computational cost and simpler wavefields (scalar wavefields). However, conventional finite-difference solutions of the fourth order pseudo-acoustic wave equation for transversely isotropic (TI) media suffer from unwanted shear wave artefacts (SV waves), and stability limitations for negative ellipticity (η < 0). Alkhalifah (2000) provided a simple solution to remove the unwanted SV waves by placing the sources into isotropic areas. Since then many other methods including generalized-screen propagators (Han & Wu 2003), the auxiliary functions (Zhou et al.2006), the decoupled method (Liu et al.2009), the finite-difference frequency-domain method (Operto et al.2009) and recently proposed low-rank approximations (Fomel et al.2013; Song et al.2013; Sun et al.2015) were introduced to solve the pseudo-acoustic wave equations. Among all these methods, the low-rank approximation approach provides a good balance between the computational cost and accuracy requirements. This pseudo-analytical method optimally selects reference velocities and weights necessary to describe the wave propagation in the Fourier domain. However, the more complex the model (including strong inhomogeneity and anisotropy), the more multidimensional fast Fourier transforms (FFTs) are required per time step. In this case, the approach becomes relatively expensive. Wu & Alkhalifah (2014) proposed an optimized expansion method to reduce the required FFTs per time step. Still, the dependency on the complexity of the model limits their practical applications in complex anisotropic media. A novel way to reduce the dependency of the computational cost on model anisotropy complexity is to decompose the total extrapolated wavefield into an isotropic background part plus a correction term (Alkhalifah et al.2013; Xu & Zhou 2014; Zhang & Alkhalifah 2016). The background wavefield can be calculated using any conventional numerical method, and then the correction term is calculated using the background wavefield and anisotropy. Using an effective-isotropic model to represent the complex anisotropy for wavefield extrapolation provides a way to reduce the dependency of the computational cost on the level of anisotropy. To be specific, effective models have been handled for this purpose in two ways: (1) A static model approach: generating an equivalent velocity model by solving the anisotropic eikonal equation (Alkhalifah et al.2013; Waheed & Alkhalifah 2015). This approach solves the isotropic acoustic wave equation and is highly efficient, but only accurate for the first arrivals as the eikonal solution dictates. (2) A dynamic implementation: applying the dispersion relation for TI media to correct the isotropic wavefield in each extrapolation step by decomposing the original pseudo-differential operator into a Laplacian and a scalar operator (Xu & Zhou 2014). The Laplacian operator can be implemented using existing numerical methods, and the scalar operator can be calculated from the background wavefields with some approximations. The latter approach corrects for the kinematics of later arrivals as well, with the additional cost of computing the spatial derivatives of current wavefields. Both methods suffer from amplitude distortions. For the decomposition approach, Xu et al. (2015) proposed an elliptic decomposition (in which the background wavefield is no longer isotropic) of the dispersion equation, which provides more balanced amplitude. We propose a new spherical decomposition method and implement it using the isotropic low-rank method. The new decomposition provides more balanced amplitudes than the previous spherical decomposition. The computational complexity of our proposed method has the same order as the isotropic low-rank approximations, and it is independent of the complexity of the anisotropy. 2 THEORY The acoustic anisotropic wave equation for VTI media was derived by setting the shear wave velocities equal zero along the symmetry axis direction (Alkhalifah 1998). Alkhalifah (2000) then derived a fourth-order kinematically accurate acoustic VTI equation based on the dispersion relation   \begin{eqnarray} k_z^2 = \frac{v_{n}^2}{v_{p0}^2}\left(\frac{\omega ^2}{v_{n}^2} - \frac{\omega ^2 k_h^2}{\omega ^2-2v_n^2\eta k_h^2 }\right), \end{eqnarray} (1)where kz and kh denote the vertical and horizontal wavenumbers ($$k_h^2 = k_x^2+k_y^2$$), respectively, vp0 and $$v_n=v_{p0}\sqrt{1+2\delta }$$ are vertical P-wave and normal moveout velocities, and η is the ellipticity coefficient. The corresponding quasi-P wave equation for VTI media (Alkhalifah 2000) can be written as   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-\frac{v_{p0}^2}{2} \left((1+2\epsilon )\frac{\partial ^2 }{\partial h^2}+\frac{\partial ^2 }{\partial z^2}+Q \right)u = 0, \end{eqnarray} (2)where vp0 is vertical P-wave velocity, ε and δ are Thomsen's parameters. The differential operator $$\frac{\partial ^2}{\partial h^2}$$ is defined by $$\frac{\partial ^2}{\partial h^2} =\frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$$, where {x, y, z} define the location in Cartesian representation of space coordinates. In addition, we have a pseudo-differential term given by   \begin{eqnarray} Q=\sqrt{\left((1+2\epsilon )\frac{\partial ^2}{\partial h^2}+\frac{\partial ^2}{\partial z^2}\right)^2-8(\epsilon -\delta )\frac{\partial ^2}{\partial h^2}\frac{\partial ^2}{\partial z^2}}. \end{eqnarray} (3) The general principle behind the effective-model method is to decompose the original pseudo-differential equation into an isotropic-background equation and a corresponding correction term. The background equation can be solved efficiently by existing numerical methods (i.e. finite-difference); the correction term has an analytical expression, which can be calculated directly. Xu & Zhou (2014) proposed a spherical decomposition of eq. (2),   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-v_{p0}^2\nabla ^2 u - v_{p0}^2 \nabla \cdot (\triangle S \nabla u) = 0, \end{eqnarray} (4)where they define the ▵S operator as   \begin{eqnarray} \triangle S=\frac{1}{2}\left(\sqrt{\left((1+2\epsilon )n_{h}^2+n_{z}^2 \right)^2-8(\epsilon -\delta )n_{h}^2 n_z^2}+2\epsilon n_{h}^2 -1 \right). \end{eqnarray} (5)Here, vector $$\vec{n}=(n_x,n_y,n_z)$$ is the unit vector of phase direction ($$n_h^2=n_x^2+n_y^2$$), and can be approximated using the gradients of wavefields:   \begin{eqnarray} \vec{n}=\frac{\vec{k}}{|\vec{k}|} \approx \frac{\nabla u}{|\nabla u|}, \end{eqnarray} (6)where $$\vec{k}$$ is the vector of spatial wavenumbers and u is the current wavefield. For this approach, the first two terms on the left-hand side of eq. (4) define the background wave equation, which can be computed using regular methods. The third term is the correction term which incorporates anisotropic effects into the isotropic background wavefields. However, the correction ability of ▵S is limited, and it generates unbalanced amplitudes when the background wavefield is far from the true one. To overcome the problem, we define a more accurate background wave equation with an effective velocity of ve in all directions and derive the corresponding correction term,   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-v_{e}^2\nabla ^2 u - v_{e}^2 \triangle S \nabla ^2 u = 0. \end{eqnarray} (7)Assuming the model is homogeneous, the time-marching scheme is given by (Etgen & Brandsberg-Dahl 2009; Fomel et al.2013),   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) \int { \widehat{u} ({\bf k}, t) e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]} \text{d}{\bf k}}, \end{eqnarray} (8)with   \begin{eqnarray} \triangle S = \frac{\sqrt{\left((1+2\epsilon )n_{h}^2+n_{z}^2 \right)^2-8(\epsilon -\delta )n_{h}^2 n_z^2}-2\epsilon n_{z}^2 -1}{2(1+\epsilon )}. \end{eqnarray} (9) The effective background velocity $$v_e=v_{p0}\sqrt{1+\epsilon }$$ is averaged velocity ($$v_{p0} < v_{p0}\sqrt{1+\epsilon } < v_{p0}(1+\epsilon )$$) of the anisotropic medium. Eq. (7) is equivalent to the original quasi-P wave equation as well as the previous spherical decomposition shown in eqs (2) and (4), respectively. The use of ve guarantees that the background term in eq. (7) is closer to the exact quasi-P operator, and thus, the required corrections are overall smaller and more immune from the directional amplitude errors in eq. (6). We use an isotropic low-rank approximation to solve the background part of eq. (7). The extrapolation term $$e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]}$$ in eq. (8) can be written in a matrix form (Fomel et al.2013),   \begin{eqnarray} W({\bf x}, {\bf k}) = e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]}. \end{eqnarray} (10)For a given ▵t, the propagation matrix in eq. (10) can be approximated by a separated representation,   \begin{eqnarray} W({\bf x}, {\bf k}) \approx \sum _{m=1}^{M} \sum _{n=1}^{N} W({\bf x}, {\bf k}_m) a_{mn} W({\bf x}_n, {\bf k}). \end{eqnarray} (11)Here M and N denote the number of selected representative wavenumbers and spatial points in the model, respectively. Eq. (8) in a low-rank approximation form is given by   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) \sum _{m=1}^{M} W({\bf x}, {\bf k}_m) \left( \sum _{n=1}^{N} a_{mn} \left( \int \widehat{u} ({\bf k}, t) e^{{\bf k} \cdot {\bf x}} W({\bf x}_n, {\bf k} ) \text{d}{\bf k} \right) \right). \end{eqnarray} (12)The integral over k can be performed using N inverse FFTs per time step. With a fixed approximation error, the number N depends on the complexity of the model. For a complex isotropic model, usually 2 or 3 FFTs (N = 2 or 3) are usually required per time step, however, for an inhomogeneous anisotropic model, more than 7 or 8 FFTs are required for the same error tolerance. In some extreme cases, the number can go up to 18 (Casasanta et al.2015). Our proposed method keeps the background model isotropic, and thus, only requires 2 or 3 FFTs for low-rank extrapolation per time step even for anisotropic media. However, to calculate the gradient of background wavefield, one additional pair of FFTs is needed. The total computational cost is still much less than that required for solving the anisotropic equation directly using the regular anisotropic low-rank method. An exponential damping absorbing boundary condition is used to simulate the wavefields in a boundless space (Kosloff & Kosloff 1986). In our application, the extrapolation equation with an absorbing boundary condition is given by   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) D({\bf x})\sum _{m=1}^{M} W({\bf x}, {\bf k}_m) \left( \sum _{n=1}^{N} a_{mn} \left( \int \widehat{u} ({\bf k}, t) e^{{\bf k} \cdot {\bf x}} W({\bf x}_n, {\bf k} )\text{d}{\bf k} \right) \right). \end{eqnarray} (13)Here D(x) = e−αx decays the wavefield in the predefined absorbing areas. In practice, the damping operator is applied to the previous and current wavefields per time step. 2.1 Analysis of the correction operator The correction factor, ▵S, compensates for the inadequacy of the background wavefield in representing the anisotropic model. It depends on the anisotropy and the wave propagation directions approximated by the wavefield gradients (eq. 6). Xu et al. (2015) showed the amplitude inadequacy of this correction, so a more accurate background wavefield was needed to calculate ▵S. Our proposed method provides a more accurate background wavefield than the original one and can be implemented using the isotropic low-rank approach. Fig. 1 illustrates the correction process of the ▵S operator with the help of eikonal solutions. The red line indicates the exact dispersion curve (wave front) of a homogeneous VTI model (with parameters given by the first numerical test). The green curve depicts the background wavefield of the original spherical decomposition and the blue one depicts the background wavefield of the proposed method. The background wavefield of the proposed method is overall closer to the exact one. Different from the previous corrections, the proposed corrections are milder and thus introduce fewer distortions in amplitudes. Fig. 2 plots ▵S in eq. (7), and as expected, it is divided into four parts: two in the vertical direction and two in the horizontal direction with contrary polarities. The polarities in Fig. 2 denote the directions (faster or slower) of the corrections and the magnitudes denote the extent of corrections. The maximum magnitudes of the corrections are almost the same (∼0.002 in absolute value) in the vertical and horizontal directions since the background wavefields have nearly the same distance from the exact wavefields. In practice, wavefields can be zeros and have multiple branches at certain points, which can introduce instabilities and distortions in the correction term. To overcome these problems, an expensive wavefield decomposition into angular components at each time step or a strong regularization can be used. As an efficient but approximate alternative, we smoothed the calculated gradients with a small windowed Gaussian filter. Other strategies such as shaping regularization which regards the division as an inverse problem are also applicable (Fomel 2007). Figure 1. View largeDownload slide Comparison of dispersion curves: the exact one (red), the original background (green) and the proposed background (blue). The traveltime of the original curve is far from the exact in the horizontal direction and thus has large distortion in amplitudes. The proposed is overall closer to the exact. Figure 1. View largeDownload slide Comparison of dispersion curves: the exact one (red), the original background (green) and the proposed background (blue). The traveltime of the original curve is far from the exact in the horizontal direction and thus has large distortion in amplitudes. The proposed is overall closer to the exact. Figure 2. View largeDownload slide The correction term of the proposed method. It is divided into four parts with different polarities as expected. Polarities and magnitudes govern the direction (faster or slower) and the distance of corrections. Figure 2. View largeDownload slide The correction term of the proposed method. It is divided into four parts with different polarities as expected. Polarities and magnitudes govern the direction (faster or slower) and the distance of corrections. In summary, for each time step, the proposed method takes the following steps: Solve the first two terms of eq. (7) by an isotropic low-rank approximation. Calculate the local direction of phase velocity from current wavefields (eq. 6). Calculate the correction operator ▵S by eq. (9). Apply the correction operator to update current wavefields. 3 NUMERICAL EXAMPLES In the following tests, we extrapolate the background wavefield using the isotropic low-rank approximation. The correction term (eq. 9) is calculated from the anisotropy and the derivatives of current isotropic wavefield. Instead of using the regular finite-difference method, we use the Fourier method to calculate the derivatives and thus one additional pair of FFTs is needed in the proposed method. 3.1 Homogeneous VTI model We start with a simple homogeneous VTI example. The vertical P-wave velocity is 1500 m s−1 and the Thomsen's parameters are ε = 0.25 and δ = 0.1. The spatial and time samplings of the wavefield are 6.096 m and 1 ms, respectively. The injected source is a Ricker wavelet with a peak frequency of 30 Hz and is located at the centre of the model. Figs 3(a)–(c) show snapshots of the wavefields at 0.5 s, which are computed by previous spherical decomposition (eq. 4), the proposed spherical decomposition (eq. 7) and the anisotropic low-rank method, respectively. Here we take the anisotropic low-rank method as the true solution and compare it with the two decomposition methods. The snapshot (Fig. 3a) shows that waveforms calculated from the previous spherical decomposition method are distorted in the horizontal direction. The proposed method (Fig. 3b) suffers less distortion and is closer to the solution of anisotropic low-rank approach (Fig. 3c). To compare the waveforms more clearly, we extract two slices: the vertical and horizontal ones from Fig. 3. All the solutions have the same traveltime. However, the previous spherical decomposition has relatively larger and smaller amplitudes in the vertical and horizontal directions respectively. Waveforms of the proposed decomposition have more balanced amplitudes in different directions (vertical and horizontal) as shown in Fig. 4. We point out that, for a homogeneous medium, the decomposition method doesn't reduce the computational cost compared with the anisotropic low-rank method. Figure 3. View largeDownload slide Snapshots of the wavefield at 0.5 s. (a) The original spherical decomposition method, (b) our proposed spherical decomposition method and (c) the anisotropic low-rank method. The proposed method provides more balanced amplitude in the horizontal direction. Figure 3. View largeDownload slide Snapshots of the wavefield at 0.5 s. (a) The original spherical decomposition method, (b) our proposed spherical decomposition method and (c) the anisotropic low-rank method. The proposed method provides more balanced amplitude in the horizontal direction. Figure 4. View largeDownload slide Waveform slices in Fig. 3. (a) In depth direction and (b) in horizontal direction. The proposed method performs better in both vertical and horizontal directions compared with the previous decomposition method. Figure 4. View largeDownload slide Waveform slices in Fig. 3. (a) In depth direction and (b) in horizontal direction. The proposed method performs better in both vertical and horizontal directions compared with the previous decomposition method. 3.2 Modified Hess VTI model We next verify the effectiveness of our proposed method using the modified Hess VTI model shown in Fig. 5. A Ricker wavelet with a peak frequency of 30 Hz is injected as an explosive source. Space and time sampling of the wavefield are 6.096 m and 1 ms respectively. To verify the accuracy of the proposed method, we first deploy the source at the centre of the model and then compare the traveltimes (red lines) with the solution from the anisotropic eikonal equation (Fomel 2004). The red line in Fig. 6 is calculated by solving the VTI eikonal equation, which is regarded as the true solution here. The original spherical decomposition in Fig. 6(a) fails to generate kinematically accurate wave fronts in complex areas (as indicated by the arrow). The wave front of the proposed decomposition method in Fig. 6 fits the true solution better in the same area. This example shows that the proposed method can not only provide more balanced amplitudes (dynamic information) but also more accurate traveltimes (kinematic information). As discussed in the theory part, that is because the background wavefield of the proposed method is closer to the true solution. For the same order of decomposition errors (∼10−6 in this example), the proposed method requires three inverse FFTs plus one pair of FFTs per time step. However, the regular anisotropic low-rank decomposition needs 8 inverse FFTs per time step (one forward FFT is required by both of them). For more advanced anisotropic models (such as orthorhombic media), the difference will be even larger since the proposed method only needs an isotropic low-rank extrapolation regardless of the complexity of the anisotropy. Figure 5. View largeDownload slide Modified Hess VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. The dipping faults are quite common in the real subsurface and its existing increases the difficulty of imaging. Figure 5. View largeDownload slide Modified Hess VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. The dipping faults are quite common in the real subsurface and its existing increases the difficulty of imaging. Figure 6. View largeDownload slide Comparison of traveltimes. Red lines indicate the solution of the anisotropic eikonal equation. (a) Original spherical decomposition, (b) the proposed method. The arrow indicates the position where the original decomposition fails but the proposed method performs better. Figure 6. View largeDownload slide Comparison of traveltimes. Red lines indicate the solution of the anisotropic eikonal equation. (a) Original spherical decomposition, (b) the proposed method. The arrow indicates the position where the original decomposition fails but the proposed method performs better. We next apply the proposed method in a Reverse Time Migration (RTM) algorithm (Etgen et al.2009; Liu et al.2011). We use the same wavelet and sampling as the previous test, and we evenly deploy 70 shots and 150 receivers on the surface of the model. The reflectivity model (Fig. 7a) is regarded as the true solution. The existing dipping layers in the model can generate reflections with different directions of propagation, which can be used to determine whether the proposed method can provide balanced amplitudes versus propagation angles. Fig. 7(b) shows the migration image using the proposed extrapolation method. One possible weakness that limits the application of the effective-model method is that unbalanced amplitudes can introduce errors in the migration images and thus lead to misinterpretations of the subsurface. The example shows that the proposed method can image structures correctly in depth with balanced amplitudes compared to the true reflectivity in Fig. 7(a). Figure 7. View largeDownload slide Migration image. (a) Reflectivity model derived from the velocity model and (b) migration image using the proposed extrapolation method. The migration image is close to the true reflectivity. Figure 7. View largeDownload slide Migration image. (a) Reflectivity model derived from the velocity model and (b) migration image using the proposed extrapolation method. The migration image is close to the true reflectivity. 3.3 Modified BP VTI Model To further verify the effectiveness of the proposed method, we conduct an anisotropic RTM on the modified BP VTI model (set the anisotropy angles equal to zero). An isotropic RTM is attached as a reference. The observed data are generated using the original anisotropic low-rank method (solving eq. 2). We use the same geometry as the second example except that the time and spatial sampling of the wavefield are 1 ms and 6.25 m, respectively. Fig. 8 shows the vertical P-wave velocity and the corresponding anisotropy. In this example, the shape of the structures in the anisotropic parameter models is different from the P-wave velocity shape. The velocity model (Fig. 8a) is smoother than the anisotropic models (Figs 8b and c), and also has fewer interfaces than the anisotropic models. We build a reflector model from the derivatives of the δ model, shown in Fig. 9(a). The anisotropic RTM of the proposed method is shown in Fig. 9(b). Aside from the boundary artefacts, all the reflectors have been imaged correctly. However, reflectors of the isotropic migration image are not focused or located at the correct depth, as the arrows indicate. Because anisotropy in the real earth is very common, an anisotropic migration is necessary to image the reflectors correctly. Figure 8. View largeDownload slide Modified BP VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. Note that the structures of the anisotropy are different from the ones of P-wave velocity. Figure 8. View largeDownload slide Modified BP VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. Note that the structures of the anisotropy are different from the ones of P-wave velocity. Figure 9. View largeDownload slide Reflectors of Modified BP VTI model. (a) Reflectivity derived from the model of δ, (b) the anisotropic migration image of the proposed method and (c) the isotropic migration image. The anisotropic migration using the proposed modelling method improves the quality of image. Figure 9. View largeDownload slide Reflectors of Modified BP VTI model. (a) Reflectivity derived from the model of δ, (b) the anisotropic migration image of the proposed method and (c) the isotropic migration image. The anisotropic migration using the proposed modelling method improves the quality of image. 4 DISCUSSION Computational efficiency and accuracy always compete with each other, especially in wavefield extrapolations. The regular anisotropic low-rank approximation can simulate kinematically accurate wavefields with a relatively high computational cost, and the cost will dramatically increase for large 3-D orthorhombic media. The effective-model method decomposes the complex wavefields in anisotropic media into an isotropic background and an anisotropic correction part, which improves the computational efficiency with a loss of some accuracy. This decomposition method was considered to be kinematically accurate; however, our second example shows that the method can also have large errors in traveltimes for complex models. We find that the errors depend on the difference between the background wavefield and the true wavefield (the perturbation). To overcome this problem, we propose an ‘averaged’ velocity model to calculate the background wavefield. Xu et al. (2015) proposed to utilize two different velocities in horizontal and vertical directions, which can be taken as a simplification of a higher level anisotropy. This decomposition can provide more accurate extrapolation results with an increase of computational cost (especially for tilted media) since the background wavefield is no longer isotropic. The number of FFTs required by the isotropic low-rank approximation proposed here can be reduced further using the optimized expansion method introduced by Wu & Alkhalifah (2014). The proposed method has an approximation (eq. 6) when calculating the phase directions. At the locations which have many wavefield branches, wavefield gradients are not well defined. To be successful, the proposed method will need either (1) an expensive wavefield decomposition into angular components at each time step, or (2) a great deal of regularization (as Poynting vectors require, with attendant accuracy). Considering the computation efficiency, we chose the second approach in practice. The proposed method can be extended to more advanced anisotropic media such as orthorhombic media (see Appendix A). Since its computational cost is not dependent on the complexity of anisotropy, the proposed method can be more efficient in that case. The acoustic anisotropic approximation itself can only simulate kinematically accurate P waves, so it also has limitations in applications such as full waveform inversion (FWI). Current waveform inversion algorithms in practice are mainly utilizing kinematic information (traveltimes or phases), which means the proposed method still has a potential to be the modelling engine of FWI (Plessix & Cao 2011; Wu & Alkhalifah 2016; Wang & Tsvankin 2016; Li et al.2017; Zhang & Alkhalifah 2017). This is because, considering the complexity of the real earth, exactly fitting the dynamic information (waveforms) is usually not practical for field data sets. The proposed method is applicable for the cases where P-wave energy is dominant. For such data sets (e.g. marine data), the acoustic anisotropic approximation is capable of simulating wave propagation in the real earth. For data sets with strong converted waves, modified objective functions such as cross-correlation based ones can be used to avoid over fitting the converted waves, which requires both P-wave and S-wave propagators. 5 CONCLUSIONS We propose a new spherical decomposition approach to represent the effective model for wavefield extrapolation in anisotropic media using isotropic low-rank extrapolation. The proposed method is free of SV-wave artefacts and provides the correct phase representation with more balanced amplitudes. The success of the proposed method is attributed to the fact that the proposed background wavefield is overall closer to the true solution than the previous spherical decomposition. The background isotropic model is a kind of ‘averaged’ one between the isotropic and anisotropic models, and thus, can generate more accurate background wavefields. The most attractive feature of this approach is the independence of anisotropy on its computational cost. The only additional computational cost added to that of the isotropic implementation is the additional pairwise Fourier transforms required by the wavefield gradient calculations. The proposed method has the same limitation as of the acoustic anisotropic approximation, which is that we can only generate P-wave fields with accurate traveltimes. Our proposed method can be used as a modelling engine in the phase- or traveltime-objected inverse problems. ACKNOWLEDGEMENTS We thank Nabil Masmoudi and Hui Wang for their helpful discussions, and Rene-Edouard Plessix, Samuel Gray and one anonymous reviewer for helpful reviews. We thank KAUST for its support and specifically the seismic wave analysis group members for their valuable insights. The published results are reproducible from the open source software Madagascar (www.ahay.org). The research was partly funded by the National Natural Science Foundation of China (Grant No. 41730425). REFERENCES Alkhalifah T., 1998. Acoustic approximations for processing in transversely isotropic media, Geophysics  63( 2), 623– 631. Google Scholar CrossRef Search ADS   Alkhalifah T., 2000. An acoustic wave equation for anisotropic media, Geophysics  65( 4), 1239– 1250. Google Scholar CrossRef Search ADS   Alkhalifah T., 2003. An acoustic wave equation for orthorhombic anisotropy, Geophysics  68( 4), 1169– 1172. Google Scholar CrossRef Search ADS   Alkhalifah T., Ma X., bin Waheed U., Zuberi M., 2013. Efficient anisotropic wavefield extrapolation using effective isotropic models, in 75th EAGE Conference & Exhibition Incorporating SPE EUROPEC 2013  London, UK. Casasanta L., Xue Z., Gray S., 2015. Converted wave RTM using lowrank wavefield extrapolation, in 77th EAGE Conference and Exhibition 2015  Madrid, Spain. Etgen J., Gray S.H., Zhang Y., 2009. An overview of depth imaging in exploration geophysics, Geophysics  74( 6), WCA5– WCA17. Google Scholar CrossRef Search ADS   Etgen J.T., Brandsberg-Dahl S., 2009. The pseudo-analytical method: application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation, in SEG Technical Program Expanded Abstracts 2009  pp. 2552– 2556, Society of Exploration Geophysicists. Fomel S., 2004. On anelliptic approximations for qP velocities in VTI media, Geophys. Prospect.  52( 3), 247– 259. Google Scholar CrossRef Search ADS   Fomel S., 2007. Shaping regularization in geophysical-estimation problems, Geophysics  72( 2), R29– R36. Google Scholar CrossRef Search ADS   Fomel S., Ying L., Song X., 2013. Seismic wave extrapolation using lowrank symbol approximation, Geophys. Prospect.  61( 3), 526– 536. Google Scholar CrossRef Search ADS   Han Q., Wu R.-S., 2003. One-way dual-domain propagators for scalar P-wave in VTI media, in SEG Technical Program Expanded Abstracts 2003  pp. 157– 160, Society of Exploration Geophysicists. Kosloff R., Kosloff D., 1986. Absorbing boundaries for wave propagation problems, J. Comput. Phys.  63( 2), 363– 376. Google Scholar CrossRef Search ADS   Li V., Wang H., Tsvankin I., Díaz E., Alkhalifah T., 2017. Inversion gradients for acoustic VTI wavefield tomography, Geophysics  82( 4), 1– 72. Liu F., Morton S.A., Jiang S., Ni L., Leveille J.P., 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media, in SEG Technical Program Expanded Abstracts 2009  pp. 2844– 2848, Society of Exploration Geophysicists. Liu Y., Chang X., Jin D., He R., Sun H., Zheng Y., 2011. Reverse time migration of multiples for subsalt imaging, Geophysics  76( 5), WB209– WB216. Google Scholar CrossRef Search ADS   Operto S., Virieux J., Ribodetti A., Anderson J.E., 2009. Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic (TTI) media, Geophysics  74( 5), T75– T95. Google Scholar CrossRef Search ADS   Plessix R.-E., Cao Q., 2011. A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium, Geophys. J. Int.  185( 1), 539– 556. Google Scholar CrossRef Search ADS   Song X., Alkhalifah T., 2013. Modeling of pseudoacoustic P -waves in orthorhombic media with a low-rank approximation, Geophysics, 78( 4), C33– C40. Song X., Fomel S., Ying L., 2013. Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation, Geophys. J. Int.  193 960– 969. Google Scholar CrossRef Search ADS   Sun J., Fomel S., Ying L., 2015. Low-rank one-step wave extrapolation for reverse time migration, Geophysics  81( 1), S39– S54. Google Scholar CrossRef Search ADS   Waheed U.b., Alkhalifah T., 2015. Vertical elliptic operator for efficient wave propagation in TTI media, in SEG Technical Program Expanded Abstracts 2015  pp. 3560– 3564, Society of Exploration Geophysicists. Wang H., Tsvankin I., 2016. Feasibility of waveform inversion in acoustic orthorhombic media, in SEG Technical Program Expanded Abstracts 2016  pp. 311– 316, Society of Exploration Geophysicists. Wu Z., Alkhalifah T., 2014. The optimized expansion based low-rank method for wavefield extrapolation, Geophysics  79( 2), T51– T60. Google Scholar CrossRef Search ADS   Wu Z., Alkhalifah T., 2016. Waveform inversion for acoustic VTI media in frequency domain, in SEG Technical Program Expanded Abstracts 2016  pp. 1184– 1189, Society of Exploration Geophysicists. Xu S., Zhou H., 2014. Accurate simulations of pure quasi-P-waves in complex anisotropic media, Geophysics  79( 6), T341– T348. Google Scholar CrossRef Search ADS   Xu S., Tang B., Mu J., Zhou H., 2015. Elliptic decomposition of quasi-P wave equation, in 77th EAGE Conference and Exhibition 2015  Madrid, Spain. Zhang Z., Alkhalifah T., 2016. Efficient quasi-P wavefield extrapolation using an isotropic lowrank approximation, in 78th EAGE Conference and Exhibition 2016  Vienna, Austria. Zhang Z.-D., Alkhalifah T., 2017. Full waveform inversion using oriented time-domain imaging method for vertical transverse isotropic media, Geophys. Prospect.  doi:10.1111/1365-2478.12560. Zhou H., Zhang G., Bloor R., 2006. An anisotropic acoustic wave equation for VTI media, in 68th EAGE Conference and Exhibition incorporating SPE EUROPEC 2006  Vienna, Austria. APPENDIX A: The acoustic approximation of setting Shear-wave velocities to zero is also valid for orthorhombic media (Alkhalifah 2003). Song & Alkhalifah (2013) successfully applied the low-rank approximation method to solve the pseudo-acoustic equation for orthorhombic media. Our proposed method is also applicable to orthorhombic media. We follow Song's derivations and extend them to the decomposition approach. By setting the velocities of shear waves (vs1, vs2 and vs0) to zero, the Christoffel equation has the following form:   \begin{eqnarray} \frac{1}{\phi ^2} {\left[\begin{array}{c@{\quad}c@{\quad}c}k_x^2v_1^2(1+2\eta _1) - 1 & \sqrt{1+2\delta }k_x k_y v_1^2(1+2\eta _1) & k_x k_z v_1 v_{p0} \\ \sqrt{1+2\delta }k_x k_y v_1^2 (1+2\eta _1) & k_y^2 v_2^2 (1+2\eta _2)-1 & k_y k_z v_2 v_{p0} \\ k_x k_z v_1 v_{p0} & k_y k_z v_2 v_{p0} & k_z^2 v_{p0}^2 -1, \end{array}\right]}, \end{eqnarray} (A1)where kx, ky and kz are wavenumbers in three directions. ϕ is the phase operator for acoustic orthorhombic media. The anisotropy are defined as $$v_{p0} = \sqrt{\frac{c_{33}}{\rho }}$$, $$v_1 = \sqrt{\frac{c_{13}(c_{13}+2c_{55})+c_{33}c_{55}}{\rho (c_{33}-c{55})}}$$, $$v_2 = \sqrt{\frac{c_{23}(c_{23}+2c_{44})+c_{33}c_{44}}{\rho (c_{33}-c{44})}}$$, $$\eta _1 = \frac{c_{11}(c_{33}-c_{55})}{2c_{13}(c_{13}+2c_{55})+2c_{33}c_{55}} - \frac{1}{2}$$, $$\eta _2 = \frac{c_{22}(c_{33}-c_{44})}{2c_{23}(c_{23}+2c_{44})+2c_{33}c_{44}} - \frac{1}{2}$$, $$\delta = \frac{(c_{12}+c_{66})^2-(c_{11}-c_{66})^2}{2c_{11}(c_{11}-c_{66})}$$. The dispersion relationship (without the frequency symbol) is given by setting the determinate of matrix A1 equal to zero:   \begin{eqnarray} &&{ -\phi ^6 + \phi ^4\big(2v_1^2\eta _1k_x^2 + v_1^2 k_x^2 +2v_2^2 \eta _2 k_y^2 + v_2^2 k_y^2 +v_{p0}^2 k_z^2 \big)+ \phi ^2(v_1^4(1+\delta )(1+2\eta _1)^2k_x^2 k_y^2 - v_2^2 v_1^2 (1+2\eta _1)(1+2\eta _2)k_x^2 k_y^2 - 2v_{p0}^2v_1^2\eta _1k_x^2 k_z^2 }\nonumber\\ &&- 2v_2^2 v_{p0}^2 \eta _2 k_y^2 k_z^2)-v_1^4v_{p0}^2(1+2\delta )(1+2\eta _1)^2k_x^2 k_y^2 k_z^2 + 2 v_1^3 v_2 v_{p0}^2 \sqrt{}(1+2\delta )(1+2\eta _1)k_x^2 k_y^2 k_z^2 -v_1^2 v_2^2 v_{p0}^2(1-4\eta _1 \eta _2) k_x^2 k_y^2 k_z^2 = 0. \end{eqnarray} (A2) To derive the governing equation for P waves in acoustic orthorhombic media, we need to find the root of the cubic polynomial (Song & Alkhalifah 2013),   \begin{eqnarray} \phi ^2 = \frac{1}{6} |-2^{\frac{2}{3}}d - \frac{2\root 3 \of {2}(a^2+3b)}{d} + 2a|, \end{eqnarray} (A3)where $$a=2v_1^2\eta _1k_x^2 + v_1^2 k_x^2 +2v_2^2 \eta _2 k_y^2 + v_2^2 k_y^2 +v_{p0}^2 k_z^2$$, $$b=v_1^4(1+\delta )(1+2\eta _1)^2k_x^2 k_y^2 - v_2^2 v_1^2 (1+2\eta _1)(1+2\eta _2)k_x^2 k_y^2 - 2v_{p0}^2v_1^2\eta _1k_x^2 k_z^2 - 2v_2^2 v_{p0}^2 \eta _2 k_y^2 k_z^2$$, $$c\,=\,v_1^4v_{p0}^2(1\,+\,2\delta )(1\,+\,2\eta _1)^2k_x^2 k_y^2 k_z^2 \,+\, 2 v_1^3 v_2 v_{p0}^2 \sqrt{(1\,+\,2\delta )}(1\,+\,2\eta _1)k_x^2 k_y^2 k_z^2 \,-\,v_1^2 v_2^2 v_{p0}^2(1\,-\,4\eta _1 \eta _2) k_x^2 k_y^2 k_z^2$$, $$d = \root 3 \of {-2a^3+3(e-9c)-9ab}$$, $$e = \sqrt{|-3b^2(a^2+4b) + 6ac(2a^2+9b) + 81c^2|}$$. The background dispersion relation (only keep ϕ) we use is given by   \begin{eqnarray} \phi _b^2 = v_e^2 (k_x^2 + k_y^2 + k_z^2). \end{eqnarray} (A4)Here $$v_e=v_{p0}\sqrt{1+\frac{\epsilon _1 + \epsilon _2}{3}}$$ is a kind of averaged velocity; ε1 and ε2 are Thomsen's parameter in two directions. By solving the equation below, the correction term can be derived,   \begin{eqnarray} (1+\triangle S)\phi _b^2 = \phi ^2. \end{eqnarray} (A5) Thus, the correction term with respect to the background wavefields (ϕb, in equation A4) is given by,   \begin{eqnarray} \triangle S = \frac{\phi ^2}{\phi _b^2}-1. \end{eqnarray} (A6)The correction operator ▵S depends on the current wavefield and the anisotropy. The decomposed equation for orthorhombic media is given by eq. (7) and then solved using eq. (13). The extension to orthorhombic media is straightforward, specifically, we change the dispersion relation to the one for orthorhombic media and derive a new correction term based on the predefined background isotropic wavefield. We verify the derivations through a homogeneous orthorhombic model. Fig. A1 compares the snapshots of an isotropic background and the proposed correction method. We can find that the proposed method can correct the traveltime of the waves and the amplitudes of the corrected wavefield are balanced. Figure A1. View largeDownload slide Snapshots. (a) Isotropic background and (b) with corrections. Figure A1. View largeDownload slide Snapshots. (a) Isotropic background and (b) with corrections. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Efficient anisotropic quasi-P wavefield extrapolation using an isotropic low-rank approximation

Loading next page...
 
/lp/ou_press/efficient-anisotropic-quasi-p-wavefield-extrapolation-using-an-SPosbKeKo4
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx543
Publisher site
See Article on Publisher Site

Abstract

Summary The computational cost of quasi-P wave extrapolation depends on the complexity of the medium, and specifically the anisotropy. Our effective-model method splits the anisotropic dispersion relation into an isotropic background and a correction factor to handle this dependency. The correction term depends on the slope (measured using the gradient) of current wavefields and the anisotropy. As a result, the computational cost is independent of the nature of anisotropy, which makes the extrapolation efficient. A dynamic implementation of this approach decomposes the original pseudo-differential operator into a Laplacian, handled using the low-rank approximation of the spectral operator, plus an angular dependent correction factor applied in the space domain to correct for anisotropy. We analyse the role played by the correction factor and propose a new spherical decomposition of the dispersion relation. The proposed method provides accurate wavefields in phase and more balanced amplitudes than a previous spherical decomposition. Also, it is free of SV-wave artefacts. Applications to a simple homogeneous transverse isotropic medium with a vertical symmetry axis (VTI) and a modified Hess VTI model demonstrate the effectiveness of the approach. The Reverse Time Migration applied to a modified BP VTI model reveals that the anisotropic migration using the proposed modelling engine performs better than an isotropic migration. Numerical modelling, Seismic anisotropy, Wave propagation 1 INTRODUCTION Seismic anisotropy has been recognized as an important factor in seismic imaging and inversion. Because of the limitations we face in our computational ability and those we face in estimating anisotropy, the acoustic approximation (Alkhalifah 1998) is still the most efficient and practical way to incorporate anisotropy into regular P-wave processing routines. Compared to solving the elastic wave equations, solving the acoustic anisotropic wave equations has two main advantages, which are less computational cost and simpler wavefields (scalar wavefields). However, conventional finite-difference solutions of the fourth order pseudo-acoustic wave equation for transversely isotropic (TI) media suffer from unwanted shear wave artefacts (SV waves), and stability limitations for negative ellipticity (η < 0). Alkhalifah (2000) provided a simple solution to remove the unwanted SV waves by placing the sources into isotropic areas. Since then many other methods including generalized-screen propagators (Han & Wu 2003), the auxiliary functions (Zhou et al.2006), the decoupled method (Liu et al.2009), the finite-difference frequency-domain method (Operto et al.2009) and recently proposed low-rank approximations (Fomel et al.2013; Song et al.2013; Sun et al.2015) were introduced to solve the pseudo-acoustic wave equations. Among all these methods, the low-rank approximation approach provides a good balance between the computational cost and accuracy requirements. This pseudo-analytical method optimally selects reference velocities and weights necessary to describe the wave propagation in the Fourier domain. However, the more complex the model (including strong inhomogeneity and anisotropy), the more multidimensional fast Fourier transforms (FFTs) are required per time step. In this case, the approach becomes relatively expensive. Wu & Alkhalifah (2014) proposed an optimized expansion method to reduce the required FFTs per time step. Still, the dependency on the complexity of the model limits their practical applications in complex anisotropic media. A novel way to reduce the dependency of the computational cost on model anisotropy complexity is to decompose the total extrapolated wavefield into an isotropic background part plus a correction term (Alkhalifah et al.2013; Xu & Zhou 2014; Zhang & Alkhalifah 2016). The background wavefield can be calculated using any conventional numerical method, and then the correction term is calculated using the background wavefield and anisotropy. Using an effective-isotropic model to represent the complex anisotropy for wavefield extrapolation provides a way to reduce the dependency of the computational cost on the level of anisotropy. To be specific, effective models have been handled for this purpose in two ways: (1) A static model approach: generating an equivalent velocity model by solving the anisotropic eikonal equation (Alkhalifah et al.2013; Waheed & Alkhalifah 2015). This approach solves the isotropic acoustic wave equation and is highly efficient, but only accurate for the first arrivals as the eikonal solution dictates. (2) A dynamic implementation: applying the dispersion relation for TI media to correct the isotropic wavefield in each extrapolation step by decomposing the original pseudo-differential operator into a Laplacian and a scalar operator (Xu & Zhou 2014). The Laplacian operator can be implemented using existing numerical methods, and the scalar operator can be calculated from the background wavefields with some approximations. The latter approach corrects for the kinematics of later arrivals as well, with the additional cost of computing the spatial derivatives of current wavefields. Both methods suffer from amplitude distortions. For the decomposition approach, Xu et al. (2015) proposed an elliptic decomposition (in which the background wavefield is no longer isotropic) of the dispersion equation, which provides more balanced amplitude. We propose a new spherical decomposition method and implement it using the isotropic low-rank method. The new decomposition provides more balanced amplitudes than the previous spherical decomposition. The computational complexity of our proposed method has the same order as the isotropic low-rank approximations, and it is independent of the complexity of the anisotropy. 2 THEORY The acoustic anisotropic wave equation for VTI media was derived by setting the shear wave velocities equal zero along the symmetry axis direction (Alkhalifah 1998). Alkhalifah (2000) then derived a fourth-order kinematically accurate acoustic VTI equation based on the dispersion relation   \begin{eqnarray} k_z^2 = \frac{v_{n}^2}{v_{p0}^2}\left(\frac{\omega ^2}{v_{n}^2} - \frac{\omega ^2 k_h^2}{\omega ^2-2v_n^2\eta k_h^2 }\right), \end{eqnarray} (1)where kz and kh denote the vertical and horizontal wavenumbers ($$k_h^2 = k_x^2+k_y^2$$), respectively, vp0 and $$v_n=v_{p0}\sqrt{1+2\delta }$$ are vertical P-wave and normal moveout velocities, and η is the ellipticity coefficient. The corresponding quasi-P wave equation for VTI media (Alkhalifah 2000) can be written as   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-\frac{v_{p0}^2}{2} \left((1+2\epsilon )\frac{\partial ^2 }{\partial h^2}+\frac{\partial ^2 }{\partial z^2}+Q \right)u = 0, \end{eqnarray} (2)where vp0 is vertical P-wave velocity, ε and δ are Thomsen's parameters. The differential operator $$\frac{\partial ^2}{\partial h^2}$$ is defined by $$\frac{\partial ^2}{\partial h^2} =\frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$$, where {x, y, z} define the location in Cartesian representation of space coordinates. In addition, we have a pseudo-differential term given by   \begin{eqnarray} Q=\sqrt{\left((1+2\epsilon )\frac{\partial ^2}{\partial h^2}+\frac{\partial ^2}{\partial z^2}\right)^2-8(\epsilon -\delta )\frac{\partial ^2}{\partial h^2}\frac{\partial ^2}{\partial z^2}}. \end{eqnarray} (3) The general principle behind the effective-model method is to decompose the original pseudo-differential equation into an isotropic-background equation and a corresponding correction term. The background equation can be solved efficiently by existing numerical methods (i.e. finite-difference); the correction term has an analytical expression, which can be calculated directly. Xu & Zhou (2014) proposed a spherical decomposition of eq. (2),   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-v_{p0}^2\nabla ^2 u - v_{p0}^2 \nabla \cdot (\triangle S \nabla u) = 0, \end{eqnarray} (4)where they define the ▵S operator as   \begin{eqnarray} \triangle S=\frac{1}{2}\left(\sqrt{\left((1+2\epsilon )n_{h}^2+n_{z}^2 \right)^2-8(\epsilon -\delta )n_{h}^2 n_z^2}+2\epsilon n_{h}^2 -1 \right). \end{eqnarray} (5)Here, vector $$\vec{n}=(n_x,n_y,n_z)$$ is the unit vector of phase direction ($$n_h^2=n_x^2+n_y^2$$), and can be approximated using the gradients of wavefields:   \begin{eqnarray} \vec{n}=\frac{\vec{k}}{|\vec{k}|} \approx \frac{\nabla u}{|\nabla u|}, \end{eqnarray} (6)where $$\vec{k}$$ is the vector of spatial wavenumbers and u is the current wavefield. For this approach, the first two terms on the left-hand side of eq. (4) define the background wave equation, which can be computed using regular methods. The third term is the correction term which incorporates anisotropic effects into the isotropic background wavefields. However, the correction ability of ▵S is limited, and it generates unbalanced amplitudes when the background wavefield is far from the true one. To overcome the problem, we define a more accurate background wave equation with an effective velocity of ve in all directions and derive the corresponding correction term,   \begin{eqnarray} \frac{\partial ^2 u}{\partial t^2}-v_{e}^2\nabla ^2 u - v_{e}^2 \triangle S \nabla ^2 u = 0. \end{eqnarray} (7)Assuming the model is homogeneous, the time-marching scheme is given by (Etgen & Brandsberg-Dahl 2009; Fomel et al.2013),   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) \int { \widehat{u} ({\bf k}, t) e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]} \text{d}{\bf k}}, \end{eqnarray} (8)with   \begin{eqnarray} \triangle S = \frac{\sqrt{\left((1+2\epsilon )n_{h}^2+n_{z}^2 \right)^2-8(\epsilon -\delta )n_{h}^2 n_z^2}-2\epsilon n_{z}^2 -1}{2(1+\epsilon )}. \end{eqnarray} (9) The effective background velocity $$v_e=v_{p0}\sqrt{1+\epsilon }$$ is averaged velocity ($$v_{p0} < v_{p0}\sqrt{1+\epsilon } < v_{p0}(1+\epsilon )$$) of the anisotropic medium. Eq. (7) is equivalent to the original quasi-P wave equation as well as the previous spherical decomposition shown in eqs (2) and (4), respectively. The use of ve guarantees that the background term in eq. (7) is closer to the exact quasi-P operator, and thus, the required corrections are overall smaller and more immune from the directional amplitude errors in eq. (6). We use an isotropic low-rank approximation to solve the background part of eq. (7). The extrapolation term $$e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]}$$ in eq. (8) can be written in a matrix form (Fomel et al.2013),   \begin{eqnarray} W({\bf x}, {\bf k}) = e^{i[{\bf k} \cdot {\bf x} + v_{e}({\bf x}, {\bf k}) |{\bf k}| \triangle t]}. \end{eqnarray} (10)For a given ▵t, the propagation matrix in eq. (10) can be approximated by a separated representation,   \begin{eqnarray} W({\bf x}, {\bf k}) \approx \sum _{m=1}^{M} \sum _{n=1}^{N} W({\bf x}, {\bf k}_m) a_{mn} W({\bf x}_n, {\bf k}). \end{eqnarray} (11)Here M and N denote the number of selected representative wavenumbers and spatial points in the model, respectively. Eq. (8) in a low-rank approximation form is given by   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) \sum _{m=1}^{M} W({\bf x}, {\bf k}_m) \left( \sum _{n=1}^{N} a_{mn} \left( \int \widehat{u} ({\bf k}, t) e^{{\bf k} \cdot {\bf x}} W({\bf x}_n, {\bf k} ) \text{d}{\bf k} \right) \right). \end{eqnarray} (12)The integral over k can be performed using N inverse FFTs per time step. With a fixed approximation error, the number N depends on the complexity of the model. For a complex isotropic model, usually 2 or 3 FFTs (N = 2 or 3) are usually required per time step, however, for an inhomogeneous anisotropic model, more than 7 or 8 FFTs are required for the same error tolerance. In some extreme cases, the number can go up to 18 (Casasanta et al.2015). Our proposed method keeps the background model isotropic, and thus, only requires 2 or 3 FFTs for low-rank extrapolation per time step even for anisotropic media. However, to calculate the gradient of background wavefield, one additional pair of FFTs is needed. The total computational cost is still much less than that required for solving the anisotropic equation directly using the regular anisotropic low-rank method. An exponential damping absorbing boundary condition is used to simulate the wavefields in a boundless space (Kosloff & Kosloff 1986). In our application, the extrapolation equation with an absorbing boundary condition is given by   \begin{eqnarray} u({\bf x}, t+\triangle t) \approx (1+\triangle S) D({\bf x})\sum _{m=1}^{M} W({\bf x}, {\bf k}_m) \left( \sum _{n=1}^{N} a_{mn} \left( \int \widehat{u} ({\bf k}, t) e^{{\bf k} \cdot {\bf x}} W({\bf x}_n, {\bf k} )\text{d}{\bf k} \right) \right). \end{eqnarray} (13)Here D(x) = e−αx decays the wavefield in the predefined absorbing areas. In practice, the damping operator is applied to the previous and current wavefields per time step. 2.1 Analysis of the correction operator The correction factor, ▵S, compensates for the inadequacy of the background wavefield in representing the anisotropic model. It depends on the anisotropy and the wave propagation directions approximated by the wavefield gradients (eq. 6). Xu et al. (2015) showed the amplitude inadequacy of this correction, so a more accurate background wavefield was needed to calculate ▵S. Our proposed method provides a more accurate background wavefield than the original one and can be implemented using the isotropic low-rank approach. Fig. 1 illustrates the correction process of the ▵S operator with the help of eikonal solutions. The red line indicates the exact dispersion curve (wave front) of a homogeneous VTI model (with parameters given by the first numerical test). The green curve depicts the background wavefield of the original spherical decomposition and the blue one depicts the background wavefield of the proposed method. The background wavefield of the proposed method is overall closer to the exact one. Different from the previous corrections, the proposed corrections are milder and thus introduce fewer distortions in amplitudes. Fig. 2 plots ▵S in eq. (7), and as expected, it is divided into four parts: two in the vertical direction and two in the horizontal direction with contrary polarities. The polarities in Fig. 2 denote the directions (faster or slower) of the corrections and the magnitudes denote the extent of corrections. The maximum magnitudes of the corrections are almost the same (∼0.002 in absolute value) in the vertical and horizontal directions since the background wavefields have nearly the same distance from the exact wavefields. In practice, wavefields can be zeros and have multiple branches at certain points, which can introduce instabilities and distortions in the correction term. To overcome these problems, an expensive wavefield decomposition into angular components at each time step or a strong regularization can be used. As an efficient but approximate alternative, we smoothed the calculated gradients with a small windowed Gaussian filter. Other strategies such as shaping regularization which regards the division as an inverse problem are also applicable (Fomel 2007). Figure 1. View largeDownload slide Comparison of dispersion curves: the exact one (red), the original background (green) and the proposed background (blue). The traveltime of the original curve is far from the exact in the horizontal direction and thus has large distortion in amplitudes. The proposed is overall closer to the exact. Figure 1. View largeDownload slide Comparison of dispersion curves: the exact one (red), the original background (green) and the proposed background (blue). The traveltime of the original curve is far from the exact in the horizontal direction and thus has large distortion in amplitudes. The proposed is overall closer to the exact. Figure 2. View largeDownload slide The correction term of the proposed method. It is divided into four parts with different polarities as expected. Polarities and magnitudes govern the direction (faster or slower) and the distance of corrections. Figure 2. View largeDownload slide The correction term of the proposed method. It is divided into four parts with different polarities as expected. Polarities and magnitudes govern the direction (faster or slower) and the distance of corrections. In summary, for each time step, the proposed method takes the following steps: Solve the first two terms of eq. (7) by an isotropic low-rank approximation. Calculate the local direction of phase velocity from current wavefields (eq. 6). Calculate the correction operator ▵S by eq. (9). Apply the correction operator to update current wavefields. 3 NUMERICAL EXAMPLES In the following tests, we extrapolate the background wavefield using the isotropic low-rank approximation. The correction term (eq. 9) is calculated from the anisotropy and the derivatives of current isotropic wavefield. Instead of using the regular finite-difference method, we use the Fourier method to calculate the derivatives and thus one additional pair of FFTs is needed in the proposed method. 3.1 Homogeneous VTI model We start with a simple homogeneous VTI example. The vertical P-wave velocity is 1500 m s−1 and the Thomsen's parameters are ε = 0.25 and δ = 0.1. The spatial and time samplings of the wavefield are 6.096 m and 1 ms, respectively. The injected source is a Ricker wavelet with a peak frequency of 30 Hz and is located at the centre of the model. Figs 3(a)–(c) show snapshots of the wavefields at 0.5 s, which are computed by previous spherical decomposition (eq. 4), the proposed spherical decomposition (eq. 7) and the anisotropic low-rank method, respectively. Here we take the anisotropic low-rank method as the true solution and compare it with the two decomposition methods. The snapshot (Fig. 3a) shows that waveforms calculated from the previous spherical decomposition method are distorted in the horizontal direction. The proposed method (Fig. 3b) suffers less distortion and is closer to the solution of anisotropic low-rank approach (Fig. 3c). To compare the waveforms more clearly, we extract two slices: the vertical and horizontal ones from Fig. 3. All the solutions have the same traveltime. However, the previous spherical decomposition has relatively larger and smaller amplitudes in the vertical and horizontal directions respectively. Waveforms of the proposed decomposition have more balanced amplitudes in different directions (vertical and horizontal) as shown in Fig. 4. We point out that, for a homogeneous medium, the decomposition method doesn't reduce the computational cost compared with the anisotropic low-rank method. Figure 3. View largeDownload slide Snapshots of the wavefield at 0.5 s. (a) The original spherical decomposition method, (b) our proposed spherical decomposition method and (c) the anisotropic low-rank method. The proposed method provides more balanced amplitude in the horizontal direction. Figure 3. View largeDownload slide Snapshots of the wavefield at 0.5 s. (a) The original spherical decomposition method, (b) our proposed spherical decomposition method and (c) the anisotropic low-rank method. The proposed method provides more balanced amplitude in the horizontal direction. Figure 4. View largeDownload slide Waveform slices in Fig. 3. (a) In depth direction and (b) in horizontal direction. The proposed method performs better in both vertical and horizontal directions compared with the previous decomposition method. Figure 4. View largeDownload slide Waveform slices in Fig. 3. (a) In depth direction and (b) in horizontal direction. The proposed method performs better in both vertical and horizontal directions compared with the previous decomposition method. 3.2 Modified Hess VTI model We next verify the effectiveness of our proposed method using the modified Hess VTI model shown in Fig. 5. A Ricker wavelet with a peak frequency of 30 Hz is injected as an explosive source. Space and time sampling of the wavefield are 6.096 m and 1 ms respectively. To verify the accuracy of the proposed method, we first deploy the source at the centre of the model and then compare the traveltimes (red lines) with the solution from the anisotropic eikonal equation (Fomel 2004). The red line in Fig. 6 is calculated by solving the VTI eikonal equation, which is regarded as the true solution here. The original spherical decomposition in Fig. 6(a) fails to generate kinematically accurate wave fronts in complex areas (as indicated by the arrow). The wave front of the proposed decomposition method in Fig. 6 fits the true solution better in the same area. This example shows that the proposed method can not only provide more balanced amplitudes (dynamic information) but also more accurate traveltimes (kinematic information). As discussed in the theory part, that is because the background wavefield of the proposed method is closer to the true solution. For the same order of decomposition errors (∼10−6 in this example), the proposed method requires three inverse FFTs plus one pair of FFTs per time step. However, the regular anisotropic low-rank decomposition needs 8 inverse FFTs per time step (one forward FFT is required by both of them). For more advanced anisotropic models (such as orthorhombic media), the difference will be even larger since the proposed method only needs an isotropic low-rank extrapolation regardless of the complexity of the anisotropy. Figure 5. View largeDownload slide Modified Hess VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. The dipping faults are quite common in the real subsurface and its existing increases the difficulty of imaging. Figure 5. View largeDownload slide Modified Hess VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. The dipping faults are quite common in the real subsurface and its existing increases the difficulty of imaging. Figure 6. View largeDownload slide Comparison of traveltimes. Red lines indicate the solution of the anisotropic eikonal equation. (a) Original spherical decomposition, (b) the proposed method. The arrow indicates the position where the original decomposition fails but the proposed method performs better. Figure 6. View largeDownload slide Comparison of traveltimes. Red lines indicate the solution of the anisotropic eikonal equation. (a) Original spherical decomposition, (b) the proposed method. The arrow indicates the position where the original decomposition fails but the proposed method performs better. We next apply the proposed method in a Reverse Time Migration (RTM) algorithm (Etgen et al.2009; Liu et al.2011). We use the same wavelet and sampling as the previous test, and we evenly deploy 70 shots and 150 receivers on the surface of the model. The reflectivity model (Fig. 7a) is regarded as the true solution. The existing dipping layers in the model can generate reflections with different directions of propagation, which can be used to determine whether the proposed method can provide balanced amplitudes versus propagation angles. Fig. 7(b) shows the migration image using the proposed extrapolation method. One possible weakness that limits the application of the effective-model method is that unbalanced amplitudes can introduce errors in the migration images and thus lead to misinterpretations of the subsurface. The example shows that the proposed method can image structures correctly in depth with balanced amplitudes compared to the true reflectivity in Fig. 7(a). Figure 7. View largeDownload slide Migration image. (a) Reflectivity model derived from the velocity model and (b) migration image using the proposed extrapolation method. The migration image is close to the true reflectivity. Figure 7. View largeDownload slide Migration image. (a) Reflectivity model derived from the velocity model and (b) migration image using the proposed extrapolation method. The migration image is close to the true reflectivity. 3.3 Modified BP VTI Model To further verify the effectiveness of the proposed method, we conduct an anisotropic RTM on the modified BP VTI model (set the anisotropy angles equal to zero). An isotropic RTM is attached as a reference. The observed data are generated using the original anisotropic low-rank method (solving eq. 2). We use the same geometry as the second example except that the time and spatial sampling of the wavefield are 1 ms and 6.25 m, respectively. Fig. 8 shows the vertical P-wave velocity and the corresponding anisotropy. In this example, the shape of the structures in the anisotropic parameter models is different from the P-wave velocity shape. The velocity model (Fig. 8a) is smoother than the anisotropic models (Figs 8b and c), and also has fewer interfaces than the anisotropic models. We build a reflector model from the derivatives of the δ model, shown in Fig. 9(a). The anisotropic RTM of the proposed method is shown in Fig. 9(b). Aside from the boundary artefacts, all the reflectors have been imaged correctly. However, reflectors of the isotropic migration image are not focused or located at the correct depth, as the arrows indicate. Because anisotropy in the real earth is very common, an anisotropic migration is necessary to image the reflectors correctly. Figure 8. View largeDownload slide Modified BP VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. Note that the structures of the anisotropy are different from the ones of P-wave velocity. Figure 8. View largeDownload slide Modified BP VTI model. (a) Vertical P-wave velocity, (b) ε and (c) δ. Note that the structures of the anisotropy are different from the ones of P-wave velocity. Figure 9. View largeDownload slide Reflectors of Modified BP VTI model. (a) Reflectivity derived from the model of δ, (b) the anisotropic migration image of the proposed method and (c) the isotropic migration image. The anisotropic migration using the proposed modelling method improves the quality of image. Figure 9. View largeDownload slide Reflectors of Modified BP VTI model. (a) Reflectivity derived from the model of δ, (b) the anisotropic migration image of the proposed method and (c) the isotropic migration image. The anisotropic migration using the proposed modelling method improves the quality of image. 4 DISCUSSION Computational efficiency and accuracy always compete with each other, especially in wavefield extrapolations. The regular anisotropic low-rank approximation can simulate kinematically accurate wavefields with a relatively high computational cost, and the cost will dramatically increase for large 3-D orthorhombic media. The effective-model method decomposes the complex wavefields in anisotropic media into an isotropic background and an anisotropic correction part, which improves the computational efficiency with a loss of some accuracy. This decomposition method was considered to be kinematically accurate; however, our second example shows that the method can also have large errors in traveltimes for complex models. We find that the errors depend on the difference between the background wavefield and the true wavefield (the perturbation). To overcome this problem, we propose an ‘averaged’ velocity model to calculate the background wavefield. Xu et al. (2015) proposed to utilize two different velocities in horizontal and vertical directions, which can be taken as a simplification of a higher level anisotropy. This decomposition can provide more accurate extrapolation results with an increase of computational cost (especially for tilted media) since the background wavefield is no longer isotropic. The number of FFTs required by the isotropic low-rank approximation proposed here can be reduced further using the optimized expansion method introduced by Wu & Alkhalifah (2014). The proposed method has an approximation (eq. 6) when calculating the phase directions. At the locations which have many wavefield branches, wavefield gradients are not well defined. To be successful, the proposed method will need either (1) an expensive wavefield decomposition into angular components at each time step, or (2) a great deal of regularization (as Poynting vectors require, with attendant accuracy). Considering the computation efficiency, we chose the second approach in practice. The proposed method can be extended to more advanced anisotropic media such as orthorhombic media (see Appendix A). Since its computational cost is not dependent on the complexity of anisotropy, the proposed method can be more efficient in that case. The acoustic anisotropic approximation itself can only simulate kinematically accurate P waves, so it also has limitations in applications such as full waveform inversion (FWI). Current waveform inversion algorithms in practice are mainly utilizing kinematic information (traveltimes or phases), which means the proposed method still has a potential to be the modelling engine of FWI (Plessix & Cao 2011; Wu & Alkhalifah 2016; Wang & Tsvankin 2016; Li et al.2017; Zhang & Alkhalifah 2017). This is because, considering the complexity of the real earth, exactly fitting the dynamic information (waveforms) is usually not practical for field data sets. The proposed method is applicable for the cases where P-wave energy is dominant. For such data sets (e.g. marine data), the acoustic anisotropic approximation is capable of simulating wave propagation in the real earth. For data sets with strong converted waves, modified objective functions such as cross-correlation based ones can be used to avoid over fitting the converted waves, which requires both P-wave and S-wave propagators. 5 CONCLUSIONS We propose a new spherical decomposition approach to represent the effective model for wavefield extrapolation in anisotropic media using isotropic low-rank extrapolation. The proposed method is free of SV-wave artefacts and provides the correct phase representation with more balanced amplitudes. The success of the proposed method is attributed to the fact that the proposed background wavefield is overall closer to the true solution than the previous spherical decomposition. The background isotropic model is a kind of ‘averaged’ one between the isotropic and anisotropic models, and thus, can generate more accurate background wavefields. The most attractive feature of this approach is the independence of anisotropy on its computational cost. The only additional computational cost added to that of the isotropic implementation is the additional pairwise Fourier transforms required by the wavefield gradient calculations. The proposed method has the same limitation as of the acoustic anisotropic approximation, which is that we can only generate P-wave fields with accurate traveltimes. Our proposed method can be used as a modelling engine in the phase- or traveltime-objected inverse problems. ACKNOWLEDGEMENTS We thank Nabil Masmoudi and Hui Wang for their helpful discussions, and Rene-Edouard Plessix, Samuel Gray and one anonymous reviewer for helpful reviews. We thank KAUST for its support and specifically the seismic wave analysis group members for their valuable insights. The published results are reproducible from the open source software Madagascar (www.ahay.org). The research was partly funded by the National Natural Science Foundation of China (Grant No. 41730425). REFERENCES Alkhalifah T., 1998. Acoustic approximations for processing in transversely isotropic media, Geophysics  63( 2), 623– 631. Google Scholar CrossRef Search ADS   Alkhalifah T., 2000. An acoustic wave equation for anisotropic media, Geophysics  65( 4), 1239– 1250. Google Scholar CrossRef Search ADS   Alkhalifah T., 2003. An acoustic wave equation for orthorhombic anisotropy, Geophysics  68( 4), 1169– 1172. Google Scholar CrossRef Search ADS   Alkhalifah T., Ma X., bin Waheed U., Zuberi M., 2013. Efficient anisotropic wavefield extrapolation using effective isotropic models, in 75th EAGE Conference & Exhibition Incorporating SPE EUROPEC 2013  London, UK. Casasanta L., Xue Z., Gray S., 2015. Converted wave RTM using lowrank wavefield extrapolation, in 77th EAGE Conference and Exhibition 2015  Madrid, Spain. Etgen J., Gray S.H., Zhang Y., 2009. An overview of depth imaging in exploration geophysics, Geophysics  74( 6), WCA5– WCA17. Google Scholar CrossRef Search ADS   Etgen J.T., Brandsberg-Dahl S., 2009. The pseudo-analytical method: application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation, in SEG Technical Program Expanded Abstracts 2009  pp. 2552– 2556, Society of Exploration Geophysicists. Fomel S., 2004. On anelliptic approximations for qP velocities in VTI media, Geophys. Prospect.  52( 3), 247– 259. Google Scholar CrossRef Search ADS   Fomel S., 2007. Shaping regularization in geophysical-estimation problems, Geophysics  72( 2), R29– R36. Google Scholar CrossRef Search ADS   Fomel S., Ying L., Song X., 2013. Seismic wave extrapolation using lowrank symbol approximation, Geophys. Prospect.  61( 3), 526– 536. Google Scholar CrossRef Search ADS   Han Q., Wu R.-S., 2003. One-way dual-domain propagators for scalar P-wave in VTI media, in SEG Technical Program Expanded Abstracts 2003  pp. 157– 160, Society of Exploration Geophysicists. Kosloff R., Kosloff D., 1986. Absorbing boundaries for wave propagation problems, J. Comput. Phys.  63( 2), 363– 376. Google Scholar CrossRef Search ADS   Li V., Wang H., Tsvankin I., Díaz E., Alkhalifah T., 2017. Inversion gradients for acoustic VTI wavefield tomography, Geophysics  82( 4), 1– 72. Liu F., Morton S.A., Jiang S., Ni L., Leveille J.P., 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media, in SEG Technical Program Expanded Abstracts 2009  pp. 2844– 2848, Society of Exploration Geophysicists. Liu Y., Chang X., Jin D., He R., Sun H., Zheng Y., 2011. Reverse time migration of multiples for subsalt imaging, Geophysics  76( 5), WB209– WB216. Google Scholar CrossRef Search ADS   Operto S., Virieux J., Ribodetti A., Anderson J.E., 2009. Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic (TTI) media, Geophysics  74( 5), T75– T95. Google Scholar CrossRef Search ADS   Plessix R.-E., Cao Q., 2011. A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium, Geophys. J. Int.  185( 1), 539– 556. Google Scholar CrossRef Search ADS   Song X., Alkhalifah T., 2013. Modeling of pseudoacoustic P -waves in orthorhombic media with a low-rank approximation, Geophysics, 78( 4), C33– C40. Song X., Fomel S., Ying L., 2013. Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation, Geophys. J. Int.  193 960– 969. Google Scholar CrossRef Search ADS   Sun J., Fomel S., Ying L., 2015. Low-rank one-step wave extrapolation for reverse time migration, Geophysics  81( 1), S39– S54. Google Scholar CrossRef Search ADS   Waheed U.b., Alkhalifah T., 2015. Vertical elliptic operator for efficient wave propagation in TTI media, in SEG Technical Program Expanded Abstracts 2015  pp. 3560– 3564, Society of Exploration Geophysicists. Wang H., Tsvankin I., 2016. Feasibility of waveform inversion in acoustic orthorhombic media, in SEG Technical Program Expanded Abstracts 2016  pp. 311– 316, Society of Exploration Geophysicists. Wu Z., Alkhalifah T., 2014. The optimized expansion based low-rank method for wavefield extrapolation, Geophysics  79( 2), T51– T60. Google Scholar CrossRef Search ADS   Wu Z., Alkhalifah T., 2016. Waveform inversion for acoustic VTI media in frequency domain, in SEG Technical Program Expanded Abstracts 2016  pp. 1184– 1189, Society of Exploration Geophysicists. Xu S., Zhou H., 2014. Accurate simulations of pure quasi-P-waves in complex anisotropic media, Geophysics  79( 6), T341– T348. Google Scholar CrossRef Search ADS   Xu S., Tang B., Mu J., Zhou H., 2015. Elliptic decomposition of quasi-P wave equation, in 77th EAGE Conference and Exhibition 2015  Madrid, Spain. Zhang Z., Alkhalifah T., 2016. Efficient quasi-P wavefield extrapolation using an isotropic lowrank approximation, in 78th EAGE Conference and Exhibition 2016  Vienna, Austria. Zhang Z.-D., Alkhalifah T., 2017. Full waveform inversion using oriented time-domain imaging method for vertical transverse isotropic media, Geophys. Prospect.  doi:10.1111/1365-2478.12560. Zhou H., Zhang G., Bloor R., 2006. An anisotropic acoustic wave equation for VTI media, in 68th EAGE Conference and Exhibition incorporating SPE EUROPEC 2006  Vienna, Austria. APPENDIX A: The acoustic approximation of setting Shear-wave velocities to zero is also valid for orthorhombic media (Alkhalifah 2003). Song & Alkhalifah (2013) successfully applied the low-rank approximation method to solve the pseudo-acoustic equation for orthorhombic media. Our proposed method is also applicable to orthorhombic media. We follow Song's derivations and extend them to the decomposition approach. By setting the velocities of shear waves (vs1, vs2 and vs0) to zero, the Christoffel equation has the following form:   \begin{eqnarray} \frac{1}{\phi ^2} {\left[\begin{array}{c@{\quad}c@{\quad}c}k_x^2v_1^2(1+2\eta _1) - 1 & \sqrt{1+2\delta }k_x k_y v_1^2(1+2\eta _1) & k_x k_z v_1 v_{p0} \\ \sqrt{1+2\delta }k_x k_y v_1^2 (1+2\eta _1) & k_y^2 v_2^2 (1+2\eta _2)-1 & k_y k_z v_2 v_{p0} \\ k_x k_z v_1 v_{p0} & k_y k_z v_2 v_{p0} & k_z^2 v_{p0}^2 -1, \end{array}\right]}, \end{eqnarray} (A1)where kx, ky and kz are wavenumbers in three directions. ϕ is the phase operator for acoustic orthorhombic media. The anisotropy are defined as $$v_{p0} = \sqrt{\frac{c_{33}}{\rho }}$$, $$v_1 = \sqrt{\frac{c_{13}(c_{13}+2c_{55})+c_{33}c_{55}}{\rho (c_{33}-c{55})}}$$, $$v_2 = \sqrt{\frac{c_{23}(c_{23}+2c_{44})+c_{33}c_{44}}{\rho (c_{33}-c{44})}}$$, $$\eta _1 = \frac{c_{11}(c_{33}-c_{55})}{2c_{13}(c_{13}+2c_{55})+2c_{33}c_{55}} - \frac{1}{2}$$, $$\eta _2 = \frac{c_{22}(c_{33}-c_{44})}{2c_{23}(c_{23}+2c_{44})+2c_{33}c_{44}} - \frac{1}{2}$$, $$\delta = \frac{(c_{12}+c_{66})^2-(c_{11}-c_{66})^2}{2c_{11}(c_{11}-c_{66})}$$. The dispersion relationship (without the frequency symbol) is given by setting the determinate of matrix A1 equal to zero:   \begin{eqnarray} &&{ -\phi ^6 + \phi ^4\big(2v_1^2\eta _1k_x^2 + v_1^2 k_x^2 +2v_2^2 \eta _2 k_y^2 + v_2^2 k_y^2 +v_{p0}^2 k_z^2 \big)+ \phi ^2(v_1^4(1+\delta )(1+2\eta _1)^2k_x^2 k_y^2 - v_2^2 v_1^2 (1+2\eta _1)(1+2\eta _2)k_x^2 k_y^2 - 2v_{p0}^2v_1^2\eta _1k_x^2 k_z^2 }\nonumber\\ &&- 2v_2^2 v_{p0}^2 \eta _2 k_y^2 k_z^2)-v_1^4v_{p0}^2(1+2\delta )(1+2\eta _1)^2k_x^2 k_y^2 k_z^2 + 2 v_1^3 v_2 v_{p0}^2 \sqrt{}(1+2\delta )(1+2\eta _1)k_x^2 k_y^2 k_z^2 -v_1^2 v_2^2 v_{p0}^2(1-4\eta _1 \eta _2) k_x^2 k_y^2 k_z^2 = 0. \end{eqnarray} (A2) To derive the governing equation for P waves in acoustic orthorhombic media, we need to find the root of the cubic polynomial (Song & Alkhalifah 2013),   \begin{eqnarray} \phi ^2 = \frac{1}{6} |-2^{\frac{2}{3}}d - \frac{2\root 3 \of {2}(a^2+3b)}{d} + 2a|, \end{eqnarray} (A3)where $$a=2v_1^2\eta _1k_x^2 + v_1^2 k_x^2 +2v_2^2 \eta _2 k_y^2 + v_2^2 k_y^2 +v_{p0}^2 k_z^2$$, $$b=v_1^4(1+\delta )(1+2\eta _1)^2k_x^2 k_y^2 - v_2^2 v_1^2 (1+2\eta _1)(1+2\eta _2)k_x^2 k_y^2 - 2v_{p0}^2v_1^2\eta _1k_x^2 k_z^2 - 2v_2^2 v_{p0}^2 \eta _2 k_y^2 k_z^2$$, $$c\,=\,v_1^4v_{p0}^2(1\,+\,2\delta )(1\,+\,2\eta _1)^2k_x^2 k_y^2 k_z^2 \,+\, 2 v_1^3 v_2 v_{p0}^2 \sqrt{(1\,+\,2\delta )}(1\,+\,2\eta _1)k_x^2 k_y^2 k_z^2 \,-\,v_1^2 v_2^2 v_{p0}^2(1\,-\,4\eta _1 \eta _2) k_x^2 k_y^2 k_z^2$$, $$d = \root 3 \of {-2a^3+3(e-9c)-9ab}$$, $$e = \sqrt{|-3b^2(a^2+4b) + 6ac(2a^2+9b) + 81c^2|}$$. The background dispersion relation (only keep ϕ) we use is given by   \begin{eqnarray} \phi _b^2 = v_e^2 (k_x^2 + k_y^2 + k_z^2). \end{eqnarray} (A4)Here $$v_e=v_{p0}\sqrt{1+\frac{\epsilon _1 + \epsilon _2}{3}}$$ is a kind of averaged velocity; ε1 and ε2 are Thomsen's parameter in two directions. By solving the equation below, the correction term can be derived,   \begin{eqnarray} (1+\triangle S)\phi _b^2 = \phi ^2. \end{eqnarray} (A5) Thus, the correction term with respect to the background wavefields (ϕb, in equation A4) is given by,   \begin{eqnarray} \triangle S = \frac{\phi ^2}{\phi _b^2}-1. \end{eqnarray} (A6)The correction operator ▵S depends on the current wavefield and the anisotropy. The decomposed equation for orthorhombic media is given by eq. (7) and then solved using eq. (13). The extension to orthorhombic media is straightforward, specifically, we change the dispersion relation to the one for orthorhombic media and derive a new correction term based on the predefined background isotropic wavefield. We verify the derivations through a homogeneous orthorhombic model. Fig. A1 compares the snapshots of an isotropic background and the proposed correction method. We can find that the proposed method can correct the traveltime of the waves and the amplitudes of the corrected wavefield are balanced. Figure A1. View largeDownload slide Snapshots. (a) Isotropic background and (b) with corrections. Figure A1. View largeDownload slide Snapshots. (a) Isotropic background and (b) with corrections. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial