Reverse time migration (RTM) relies on accurate wave extrapolation engines to image complex subsurface structures. To construct such operators with high efficiency and numerical stability, we have developed a one-step wave extrapolation approach using complex-valued low-rank decomposition to approximate the mixed-domain space-wavenumber wave extrapolation symbol. The low-rank one-step method involves a complex-valued phase function, which is more flexible than a real-valued phase function of two-step schemes, and thus it is capable of modeling a wider variety of dispersion relations. Two novel designs of the phase function leads to the desired properties in wave extrapolation. First, for wave propagation in inhomogeneous media, including a velocity gradient term assures a more accurate phase behavior, particularly when the velocity variations are large. Second, an absorbing boundary condition, which is propagation-direction-dependent, can be incorporated into the phase function as an anisotropic attenuation term. This term allows waves to travel parallel to the boundary without absorption, thus reducing artificial reflections at wide incident angles. Using numerical experiments, we revealed the stability improvement of a one-step scheme in comparison with two-step schemes. We observed the low-rank one-step operator to be remarkably stable and capable of propagating waves using large time step sizes, even beyond the Nyquist limit. The stability property can help to minimize the computational cost of seismic modeling or RTM. The low-rank one-step wave extrapolation also handles anisotropic wave propagation accurately and efficiently. When applied to RTM in anisotropic media, the proposed method generated high-quality images.
The task of wave extrapolation is propagation of waves in depth or time, which finds applications in seismic modeling and migration. Conventionally, it is implemented by finite differences (FD). FD methods have a low-computational cost but suffer from dispersion artifacts and instabilities (Kosloff and Baysal, 1982). Instead, numerical differentiation of space coordinates in wave equations can be implemented by using the Fourier transform, which is known as the pseudospectral method (Fornberg, 1998). Pseudospectral methods have higher accuracy and can suppress numerical dispersion (Reshef et al., 1988; Virieux et al., 2011). On the other hand, they are limited to small time steps and may be subject to dispersion because of FD approximations used for the time derivatives. Recently, alternative strategies have been proposed to propagate waves by mixed-domain space-wavenumber operators (Soubaras and Zhang, 2008; Wards et al., 2008; Etgen and Brandsberg-Dahl, 2009; Liu et al., 2009; Zhang and Zhang, 2009; Du et al., 2010; Pestana and Stoffa, 2010; Chu and Stoffa, 2011; Fomel et al., 2013; Wu and Alkhalifah, 2014). Fowler (2010b) and Du et al. (2014) refer to these methods as recursive integral time extrapolation (RITE) methods.
RITE methods are designed to make time extrapolation stable and dispersion free in heterogeneous media for large time steps, even beyond the Nyquist limit (Du et al., 2014), and they therefore are particularly suitable for reverse time migration (RTM), a depth-migration method in which waves are extrapolated in time (Baysal et al., 1983; McMechan, 1983; Whitmore, 1983; Farmer et al., 2006; Fletcher et al., 2009; Fowler et al., 2010a). RTM performs imaging by solving the two-way wave equation in the most straightforward manner compared with other methods, and therefore it is capable of handling complicated waveforms, such as prismatic waves, which helps generate images free of artifacts that are caused by approximations of the physics of wave propagation in other methods (Etgen et al., 2009; Leveille et al., 2011). The efficiency of RTM thus relies on the wave-propagation engine.
Among the different RITE approaches, the low-rank wave extrapolation (Fomel et al., 2010, 2013) distinguishes itself by its high efficiency and flexible control over approximation accuracy. Low-rank approximation has been implemented under different frameworks, including low-rank FDs and low-rank Fourier FDs (Song et al., 2013; Fang et al., 2014). A low-rank algorithm decomposes the original mixed-domain wave propagation matrix into a small set of representative spatial locations and a small set of representative wavenumbers. Similar to other spectral methods, the cost of computation per time step with the low-rank method is , where is the total size of the computational grid and is a small number (the rank of the approximation) controlling the trade-off between accuracy and efficiency.
Fomel et al. (2013) implement low-rank wave extrapolation using a two-step time marching scheme, which involves a real valued wave-propagation operator. The application of a two-step scheme is constrained by its requirement of a real-valued phase function. In this paper, we propose adopting a one-step scheme, which propagates a complex-valued wavefield using a complex phase function. Following Zhang and Zhang (2009), we show that separation of forward- and backward-propagating waves can be achieved by constructing a complex wavefield, with its imaginary part being the Hilbert transform of the real part. The complex wavefield corresponds to the analytical signal (Taner et al., 1979). In practice, the proposed one-step scheme demonstrates significantly improved stability. Its ability to extrapolate waves using large time steps should help to reduce the cost of computationally intensive tasks, such as RTM and time-domain full-waveform inversion (FWI). A complex-valued phase function is capable of incorporating modified forms of dispersion relations. In particular, we show that by including a velocity gradient term into the approximation of the phase function, the accuracy of wave extrapolation can be improved in media with large velocity variations. Additionally, we propose a propagation-direction-dependent absorbing boundary condition that can be incorporated into the complex-valued phase function. This condition behaves like an anisotropic attenuation term, which attenuates waves preferentially in the direction perpendicular to the absorbing boundary, thus reducing artificial reflections at wide incident angles. Such modifications take advantage of the complex-valued phase function in the one-step scheme. The proposed method is easily extended to anisotropic wave propagation, such as tilted transversely isotropic (TTI) or orthorhombic media, without producing unwanted residual pseudo-S-wave energy (Du et al., 2014). We use numerical examples with synthetic models to test the accuracy, efficiency, and stability of the proposed low-rank one-step wave extrapolation method. Finally, we apply low-rank one-step RTM to the BP 2007 TTI synthetic data set to test its ability of generating high-quality seismic images.
Analytical solutions in constant velocity
Variable velocity and anisotropy
As detailed in Appendix A, the unconditional numerical stability of the one-step scheme (equation 17) can be proven theoretically in the continuous case. In practice, we have observed that the low-rank approximation preserves the stability property, as will be demonstrated in the next section with numerical examples.
Higher order terms from the phase function
Direction-dependent absorbing boundary conditions
The proposed direction-dependent absorbing boundary conditions could also be incorporated into two-step-based wave extrapolation methods (Fomel et al., 2013; Song et al., 2013; Fang et al., 2014) using low-rank approximation, by separately applying the damping operator on the wavefield and its time derivative. However, because two-step schemes do not allow direct incorporation of the absorbing boundary condition into the wave propagation matrix, additional Fourier transforms will be required. The one-step implementation is more computationally efficient and straightforward.
Reverse time migration imaging conditions
In this section, we use several numerical examples to demonstrate the properties of low-rank one-step wave extrapolation.
Complex-valued low-rank approximation
We first test the accuracy of low-rank approximation applied to a wave extrapolation matrix in a 1D inhomogeneous medium. Using a simple velocity profile with a sharp velocity contrast (Figure 1), and a time step size of 0.01 s, the real and imaginary parts of the wave extrapolation matrix defined by equation 25 with only the term are plotted in Figure 2a and 2b, respectively. An accuracy threshold of leads to an approximation rank . The approximation error is plotted in Figure 3a and 3b, with the maximum error corresponding to the prescribed accuracy requirement. To see that the accuracy threshold is strict enough to guarantee kinematic accuracy, we first use an exact matrix multiplication to calculate the exact wavefield from an initial condition (Figure 4a). Next, we use low-rank wave extrapolation to compute the wavefield and calculate their difference (Figure 4b). A negligible error can be observed from the difference section, indicating the high accuracy of low-rank wave extrapolation.
We use a simple two-layer velocity model similar to the one used by Du et al. (2014) to demonstrate the stability of one-step wave extrapolation using low-rank approximation. Figure 5 shows the comparison among the stability of low-rank one-step, low-rank two-step, and fourth-order FD methods. The velocity model has a sharp contrast at the depth of 3795 m: The upper layer has a velocity of , and the lower layer has a velocity of (Figure 5a). The model is discretized on a grid with a spacing of 15 m along the horizontal and vertical directions. An explosive source, with a Ricker wavelet using a peak frequency of 16 Hz (the maximum frequency is approximately 50 Hz), is injected in the center of the model. When a time step of 2 ms is used, the classic fourth-order FD method suffers from visible dispersion artifacts (Figure 5b), whereas the one- and two-step schemes produce waves free of artifacts (Figure 5c and 5d). When a time step of 3 ms is used, the one-step scheme is stable (Figure 5e), whereas the two-step scheme starts to develop artifacts near the velocity contrast (Figure 5f). The FD method is no longer stable, and therefore it is not plotted. At 10 ms, which corresponds to the Nyquist sampling rate, the one-step scheme remains stable (Figure 5g), but the two-step scheme becomes unstable, and thus it is not plotted. Using 20 ms, the one-step scheme is still stable, but it starts to develop ringing artifacts similar to those observed by Du et al. (2014) (Figure 5h). Using a time step size of 30 ms, the ringing effects aggravate; however, the operator remains stable (Figure 5i).
Improved phase accuracy
To test if a more accurate wavefield can be obtained by one-step extrapolation with equation 24, we first use a synthetic model with a smooth velocity distribution and large velocity variations (Figure 6a). We propagate a wavefield with the Ricker-wavelet source, using a time step of 3.5 ms. The model is discretized on a grid with spacing of 5 m along the horizontal and vertical directions. The source is injected in the center of the model. Figure 6b shows the reference wavefield propagated using only the term from equation 20, but with an exceedingly small time step (0.175 ms), so that the term is negligible and the wavefield can be treated as accurate. Figure 6c shows the difference between the reference wavefield and the wavefield propagated by equation 24 after normalization. Figure 6d shows the difference between the reference wavefield and the wavefield propagated without including the term in the phase function, also after normalization. After including the velocity-gradient term, the error decreases significantly, and it is caused mainly by small amplitude differences. On the other hand, the difference calculated without using the velocity-gradient term is due mostly to the shift in phase. This observation is further supported by Figure 7, which shows a comparison between traces extracted from the 2D wavefield snapshots at . The trace calculated using the velocity-gradient term aligns with the reference trace, whereas the trace calculated without including the velocity-gradient term has a noticeable shift in space relative to the reference trace. The shift increases with the degree of velocity variation.
Next, we perform a similar experiment using a more complicated Marmousi velocity model (Figure 8a). The model is smoothed and discretized on a grid with spacing of 25 m. The reference wavefield is propagated using a time step size of 1.5 ms (Figure 8b). Figure 8c and 8d demonstrates the difference between the reference wavefield and the wavefields calculated using phase functions with and without the velocity gradient term, both using a time step size of 15 ms. Figure 9 overlays the trace from the reference wavefield extracted at with corresponding traces from the wavefields propagated using the larger time step size. The comparison shows that, even with moderate velocity variations but a relatively large time step size, the velocity gradient term can have a noticeable contribution to the phase accuracy.
In the remaining examples of this paper, the velocity gradient term was not included in the phase function for simplicity.
Absorbing boundary conditions
In the next example, we incorporate the propagation-direction-dependent absorbing boundary condition in the wave extrapolation operator to attenuate waves reaching the boundary. The absorbing boundary defined by equation 26 does not attenuate any energy if the wave propagates parallel to the boundary, and it reaches maximum damping when the wave propagates normal to the boundary. This property reduces artificial reflected energy at large incident angles. Figure 10 compares the propagation-direction-dependent absorbing boundary condition with a conventional direction-independent absorbing boundary condition (tapering) using a point-source wavefield. After normalization and clipping for display, the direction-dependent absorbing boundary condition proves to be more effective at wide incident angles. In isotropic wave propagation, we observed that the direction-dependent absorbing boundary condition may increase the rank of the low-rank approximation because of introduced anisotropic attenuation of wave propagation. One possible way to decrease the cost is to implement the absorbing boundary condition separately from the wave extrapolation process, using partial FFT (Ying and Fomel, 2009). In this way, the cost of absorbing boundaries can be marginal compared with the cost of wave extrapolation.
Wave propagation in tilted transverse isotropic media
To show that the proposed method handles anisotropic media accurately, we propagate qP- and qSV-waves in a portion of the 2007 anisotropic benchmark model from BP. The model has horizontal spacing of 50 m and vertical spacing of 25 m. As shown in Figure 11, the model exhibits strong tilted transverse isotropy (TTI). We restrict the modeled area to between 34 and 60 km. Because the original model does not include S-wave velocity, we compute a vertical S-wave velocity profile as 0.3 times the original vertical P-wave velocity.
To mimic an unbounded medium, we apply the absorbing boundary condition in equation 26. First, a time step size of 4 ms for the qP-wave and 2 ms for the qSV-wave are used. A Ricker wavelet source with a peak frequency of 16 Hz is located at and . Figures 12a and 13a illustrate wavefield snapshots taken at different times overlaid on top of each other. Because of the accumulative effect of the large time step, model complexity, and extra anisotropy introduced by the absorbing boundary conditions, the low-rank approximation took for the qP-wave and for the qSV-wave for an accuracy threshold of .
To test the effect of using large for wave propagation, we choose a series of increasing , while enforcing the same accuracy requirement (). The qP-wave wavefields (Figure 12) generated using , 10, 20, and 30 ms are almost identical. Similarly, the qSV-wave wavefields (Figure 13) generated using , 4, 10, and 20 ms have very little difference. This shows that the stability and accuracy of wave extrapolation are not compromised by very large time steps, even beyond the Nyquist limit of 10 ms.
To investigate the effects of time step size () and accuracy level () of low-rank approximation on the computational cost, which is directly controlled by the rank of the approximation, we conduct a series of experiments using different values of and for qP-wave propagation. No absorbing boundary condition is applied in this test. Table 1 lists the rank required by different time step sizes to achieve different accuracy levels. We observe that the low-rank algorithm requires a higher rank to maintain the same level of accuracy when the time step size increases. On the other hand, with the same time step size, a higher rank may be needed to achieve higher accuracy. In practice, we find that is small enough to guarantee accurate wave kinematics. Because making a larger step in time requires a smaller number of steps, the optimal time step size can be selected to be the one that minimizes the total number of FFTs required to accomplish the modeling or imaging task at hand. In practice, we found that this minimization problem usually has a straightforward solution: The larger is, the less computation is required. Therefore, the optimal solution will be the largest that satisfies other constraints, such as the imaging condition when performing RTM.
Wave propagation in 3D orthorhombic media
RTM of BP 2007 TTI data set
Finally, we test low-rank one-step RTM using the BP 2007 anisotropic benchmark data set (Figure 11). The grid spacing is 18.75 m in the horizontal direction and 12.5 m in the vertical direction. The time step size is selected to be 4 ms. We use the acoustic approximation (Alkhalifah, 1998, 2000; Fomel, 2004) to calculate the qP-wave phase velocity. The sedimentary layers, salt boundaries, anticline, and fault surfaces are clearly imaged by the low-rank one-step RTM method (Figure 16).
The modeling experiments using an isotropic two-layer model and a portion of the BP 2007 TTI benchmark model show that the proposed operator is in practice remarkably stable and accurate. When the velocity contrast is not very sharp, the low-rank one-step method is capable of propagating waves free of dispersion artifacts using time steps even larger than the Nyquist sampling limit of the source wavelet. The maximum efficiency can be achieved using the largest step size possible, which in RTM applications usually corresponds to the time sampling of the imaging condition. In contrast, conventional approaches, such as high-order FD and pseudospectral methods, are confined to small step sizes due to severe stability and accuracy constraints. This is especially true when S-wave extrapolation is required by RTM because the slow S-wave velocity at shallow depths imposes a strict constraint on the time step size of conventional methods. Low-rank one-step wave extrapolation has been recently applied to converted wave imaging by Casasanta et al. (2015) and has achieved accurate results.
Our numerical experiments confirm the advantages of low-rank one-step wave extrapolation over the two-step scheme. By extrapolating an analytical wavefield with the imaginary part related to the first derivative of the real-valued wavefield, the one-step scheme is capable of using a significantly larger time step size (Du et al., 2014). The complex phase function used by the low-rank one-step method offers additional freedom in the design of the wave extrapolation operator. In media with smoothly varying velocity, a large time step may impair accuracy using the conventional formulation. By using a more accurate expression, for example, by admitting more terms from the Taylor expansion of the phase function, the low-rank one-step wave extrapolation can achieve higher accuracy. The complex-valued extrapolation operator also allows for an effective mixed-domain absorbing boundary condition, which dampens wave energy according to the phase direction, and thus avoids artificial reflections at large incident angles. The application of a complex phase function is not limited to the proposed cases. Another possible extension is seismic modeling and imaging in viscoacoustic media (Zhu and Harris, 2014), where the one-step extrapolator can incorporate a complex-valued, decoupled dispersion relation that incorporates attenuation effects (Sun et al., 2014, 2015).
We have developed low-rank wave extrapolation using a one-step scheme. The one-step low-rank method appears to be more stable than the two-step method, and it exhibits superior stability in numerical experiments. The capability of propagating waves using large time steps can help in saving costs when performing modeling, imaging, or FWI tasks. We propose two modifications to the complex phase function, which can be accurately handled by the low-rank one-step approach. First, we show that, when the velocity-gradient term is included in the approximation of the phase function, higher accuracy can be achieved. Next, we use a one-step scheme to incorporate a propagation-direction-dependent absorbing boundary condition in the wave propagation operator, which reduces artificial reflections at wide incident angles. Numerical examples using simple models demonstrate the improved accuracy and efficiency of the proposed method. In particular, we apply the low-rank one-step operator to wave extrapolation in 2D TTI media and 3D orthorhombic media to produce wavefields free of dispersion artifacts or residual pseudo-S-wave artifacts. When low-rank one-step RTM is applied on the BP 2007 TTI benchmark data set, it produces a high-quality seismic image.
We thank Z. Wu and one anonymous reviewer for their constructive comments. We would like to also thank J. Cheng, G. Fang, P. Fowler, J. Hu, S. Li, X. Song, Z. Xue, and Y. Zhang for inspiring discussions. We thank the sponsors of the Texas Consortium for Computational Seismology for financial support. The first author was additionally supported by the Statoil Fellows Program at the University of Texas at Austin. The BP 2007 benchmark data set was created by H. Shah and released by BP Exploration Operation Co., Ltd. We thank the Texas Advanced Computing Center for providing the computational resources used in this study.
PROOF OF THE STABILITY OF THE ONE-STEP WAVE EXTRAPOLATION OPERATOR
We now apply the above lemma to . When , where is the spatial dimension, we have bounded. Therefore, for sufficiently smooth , we have , for sufficiently small . Hence, and . Because and as the Fourier transform is an isometry, we have .