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 O(NNxlogNx), where Nx is the total size of the computational grid and N 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.

Let p(x,t) be the seismic wavefield at location x and time t. The acoustic constant-density wave equation can be written as
(2t2V2(x)2)p(x,t)=0,
(1)
where V(x) is the velocity and 2 denotes the Laplacian operator.

Analytical solutions in constant velocity

When V is constant, after Fourier transform in space, the wave equation takes the form
(2t2+V2|k|2)P(k,t)=0,
(2)
where k is the spatial wavenumber and P(k,t) is the spatial Fourier transform of p(x,t):
P(k,t)=1(2π)3p(x,t)eik·xdx.
(3)
The analytical solution to equation 2 can be expressed as
P=A1ei|k|Vt+A2ei|k|Vt=P1+P2,
(4)
where P1 represents the forward-propagating wavefield, i.e., positive frequencies, and P2 represents the backward-propagating wavefield, i.e., negative frequencies. The time derivative of P has the following form:
Pt=i|k|V(P1P2).
(5)
Zhang and Zhang (2009) use the Hilbert transform to define an additional function:
Q(k,t)=1ψP(k,t)t,
(6)
where Q(k,t) is the Hilbert transform of P(k,t), and ψ=V|k|. Combining equations 4–6, P1 and P2 can be expressed as
P1=12(PiQ)
(7)
and
P2=12(P+iQ).
(8)
Equation 2 can be split into a pair of first-order equations and expressed in the following matrix form:
t[PPt]=[01ψ20][PPt].
(9)
With the help of the Hilbert transform, a more symmetric expression can be achieved as
t[PiQ]=[0iψiψ0][PiQ].
(10)
We can further decompose the first matrix on the right side as follows:
[0iψiψ0]=[1111][iψ00iψ][1/21/21/21/2].
(11)
Substituting equation 11 into equation 10, and using equations 7 and 8, we arrive at
t[PiQ]=[1111][iψ00iψ][P1P2].
(12)
In RTM, only one branch of the total wavefield is needed at one time. The two parts of wave propagation decouple according to
t[P1P2]=[iψ00iψ][P1P2].
(13)
Modeling seismic wave propagation requires the source function. Letting the source function be f(x,t), wave equation 2 can be rewritten in the following form:
(2t2+ψ2)P(k,t)=f^(k,t).
(14)
Correspondingly, equation 13 becomes
t[P1P2]=[1/21/21/21/2]{[0iψiψ0][PiQ]+[0iψf^]}=[iψ00iψ][P1P2]+[i2ψf^i2ψf^].
(15)
The application of operator i/2ψ can be implemented in either time domain or Fourier domain; it can also be directly incorporated into the definition of the source functions. For example, operator 1/ψ can be regarded as (iω/|ω|)·(1/iω), which in the time domain corresponds to cascading the Hilbert transform with the first-order integration.
In constant velocity, the forward-propagating wavefield away from the source at the next time step t+Δt can be expressed as
p1(x,t+Δt)=P1(k,t)ei[k·x+V|k|Δt]dk.
(16)

Variable velocity and anisotropy

For inhomogeneous and anisotropic media, we can use the general form of equation 16 to introduce a general phase function, which depends on k and x. We rewrite equation 16 as the following mixed-domain operator (Wards et al., 2008; Fomel et al., 2013):
p1(x,t+Δt)=P1(k,t)eiϕ(x,k,Δt)dk.
(17)
Equation 17 is kinematically correct if the phase function ϕ(x,k,t) satisfies the anisotropic eikonal equation:
ϕt=±V(x,k)|ϕ|,
(18)
where V(x,k) is the phase velocity. When the velocity is constant and isotropic, the phase function reduces to
ϕ(x,k,t)=k·x±V|k|,
(19)
and it corresponds to equation 16. In the more general case, assuming small time steps, ϕ can be expanded into the Taylor series in t (Fomel et al., 2013):
ϕ(x,k,t)k·x+ϕ1(x,k)t+ϕ2(x,k)t22+,
(20)
where
ϕ1(x,k)=V(x,k)|k|
(21)
and
ϕ2(x,k)=V(x,k)V·k.
(22)
Equation 17 corresponds to the one-step method (Zhang and Zhang, 2009; Fowler, 2010b). Fomel et al. (2013) adopt instead a two-step implementation, which uses only the ϕ1 term in equation 20 to cancel out the imaginary part of the wave extrapolation operator (Etgen and Brandsberg-Dahl, 2009):
p(x,t+Δt)+p(x,tΔt)2P(k,t)eik·xcos[V(x,k)|k|Δt]dk.
(23)
It is important to realize that the cancellation of the imaginary part relies on the fact that the phase function only contains odd-ordered terms of t, with the exception of the zeroth-order term, which corresponds to the inverse Fourier transform. Any modification to the phase function that violates this requirement cannot be easily handled by the two-step method.

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

In this paper, we adopt the one-step scheme due to its superior stability and ability to handle complex-valued phase functions. Because the two-step scheme depends on a real-valued phase function, it cannot include higher order terms from the expansion in equation 20 and thus implement more accurate phase functions. By switching to a one-step scheme, we can easily incorporate the second-order term (ϕ2) in equation 20 to achieve a more accurate wave extrapolation operator. As defined in equation 22, ϕ2 involves the gradient of the velocity model, and it can become significant when either the time step size is large or the velocity model changes rapidly. Substituting the first three terms from the Taylor series (equation 20) into equation 17, the second-order operator takes the form
p1(x,t+Δt)P1(k,t)ei[k·x+V(x,k)|k|Δt+V(x,k)V·kΔt2/2]dk.
(24)
This modification helps to increase the accuracy of wave extrapolation, especially when Δt or V is large. Note that introduction of the velocity gradient term does require the velocity model to be smoothly varying, which conforms to the usual requirement of RTM. The term V·k also leads to a certain degree of anisotropy in the phase function, which in practice may increase the numerical rank of low-rank approximation.

Direction-dependent absorbing boundary conditions

We propose another innovation using the one-step scheme by incorporating propagation-direction-dependent absorbing boundary conditions into the extrapolation operator. In absorbing boundary layers, we propose to modify the operator
W(x,k)=ei[ϕ(x,k,Δt)k·x]
(25)
into the following form:
W(x,k)abc=ei[ϕ(x,k,Δt)k·x][α(xxt)·k|k|]2,
(26)
where α is the decay parameter and xt is the location of the nearest absorbing boundary. The added term e[α(xxt)·k/|k|]2 causes exponential decay of the wavefield between the absorbing boundary and the computational boundary (domain truncation), which also depends on the wave-traveling direction. When the wave travels normal to the boundary, it will have the maximum absorbing effect because xxt is in the direction of k; on the other hand, there will be no damping when the wave is traveling parallel to the boundary because the angle between xxt and k would be zero. Allowing waves to propagate within the absorbing zone can mitigate artificial reflections from the absorbing boundary, especially at grazing incident angles. The decay term should be applied on both the real wavefield and its first-order time derivative (the imaginary part of the analytical wavefield). This is analogous to the tapering technique described by Cerjan et al. (1985). The absorbing term in the phase function can be physically interpreted as an anisotropic attenuation effect.

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.

Low-rank approximation

To implement wave propagation numerically, we use the low-rank approximation method of Fomel et al. (2013) to decompose the wave extrapolation matrix 26 into the following separated representation:
W(x,k)m=1Mn=1NW(x,km)amnW(xn,k).
(27)
The difference is that now the low-rank decomposition is implemented for complex matrices or linear operators instead of real ones. The computation of p(x,t+Δt) then becomes
p(x,t+Δt)m=1MW(x,km){n=1Namn[eixkW(xn,k)P(k,t)dk]}.
(28)
The computational cost of representation 28 is effectively equivalent to applying N inverse fast Fourier transforms (FFTs) per time step. In practice, N is a small number, typically less than five for isotropic media. The value of N may grow with increasing model complexity, such as the introduction of anisotropy. Compared with a naive straightforward implementation of equation 17, the number of floating point operations per time step is reduced from O(Nx2) to O(NNxlogNx), where Nx is the total size of the spatial grid.

Reverse time migration imaging conditions

Because the one-step wave extrapolation kernel operates in the complex domain, it requires a definition of data and reflectivity with complex values. The analytical data follow the definition of the complex wavefield in equations 7 and 9. This implies that the input data need to be Hilbert transformed along the time axis and supplied as the imaginary part before the migration process, creating an analytical signal (Taner et al., 1979). We adopt the following complex-valued cross-correlation imaging condition (Claerbout, 1985):
Ic(x)=stS¯s(x,t)Rs(x,t),
(29)
where the lowercase s denotes shots and t denotes time samples. The real part of the complex image Ic(x) is extracted and used as the final image.
Extended imaging conditions (Sava and Vasconcelos, 2011), including space-shift (Rickett and Sava, 2002; Sava and Fomel, 2003) and time-shift (Sava and Fomel, 2006) imaging conditions, can provide additional information for migration velocity analysis. The complex-valued space-shift and time-shift imaging condition for low-rank one-step RTM takes the form
Ie(x,λ,τ)=stS¯s(xλ,tτ)Rs(x+λ,t+τ),
(30)
and it can be easily implemented in the time-space domain.

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 ϕ1 term are plotted in Figure 2a and 2b, respectively. An accuracy threshold of ϵ=104 leads to an approximation rank N=4. 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.

Two-layer model

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 1500m/s, and the lower layer has a velocity of 4500m/s (Figure 5a). The model is discretized on a 400×400 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 512×512 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 ϕ1 term from equation 20, but with an exceedingly small time step (0.175 ms), so that the ϕ2 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 ϕ2 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 X=2000m. 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 376×576 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 Z=8800m 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 calculate the qP- and qSV-wave phase velocities, we use the exact formulas defined by (Gassmann, 1964):
VP2(n,x)|k|2=12[(c11+c55)k^x+(c33+c55)k^z]+12[(c11c55)k^x(c33c55)k^z]2+4(c13+c55)2k^xk^z,
(31)
and
VSV2(n,x)|k|2=12[(c11+c55)k^x+(c33+c55)k^z]12[(c11c55)k^x(c33c55)k^z]2+4(c13+c55)2k^xk^z,
(32)
where c11, c13, c33, and c55 are the density-normalized components of the elastic tensor, k is the wavenumber vector, k^x and k^z are the wavenumbers evaluated in the rotated coordinate aligned with the symmetry axis:
k^x=kxcosθ+kzsinθ,k^z=kzcosθkxsinθ,
(33)
where θ is the tilt angle measured with respect to the horizontal direction. Alternatively, the qP-wave phase velocity can be calculated using accurate approximate formulas (Alkhalifah, 1998, 2000; Fomel, 2004; Sripanich and Fomel, 2015).

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 X=47.3km and Z=0.175km. 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 N=19 for the qP-wave and N=15 for the qSV-wave for an accuracy threshold of ϵ=104.

To test the effect of using large Δt for wave propagation, we choose a series of increasing Δt, while enforcing the same accuracy requirement (ϵ=104). The qP-wave wavefields (Figure 12) generated using Δt=4, 10, 20, and 30 ms are almost identical. Similarly, the qSV-wave wavefields (Figure 13) generated using Δt=2, 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 (Δt) 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 Δt and ϵ for qP-wave propagation. No absorbing boundary condition is applied in this test. Table 1 lists the rank N 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 ϵ=104 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 Δt 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 Δt is, the less computation is required. Therefore, the optimal solution will be the largest Δt that satisfies other constraints, such as the imaging condition when performing RTM.

Wave propagation in 3D orthorhombic media

To demonstrate 3D wave propagation in tilted orthorhombic media, we use the classic model from Schoenberg and Helbig (1997), which characterizes a TI medium with vertical fractures. The density-normalized orthorhombic stiffness matrix is (Schoenberg and Helbig, 1997)
[93.62.250003.69.842.40002.252.45.937500000020000001.60000002.182].
(34)
To introduce spatial heterogeneity, we apply a moderate perturbation to the stiffness coefficients that is a function of x, as demonstrated by Figure 14 for the case of C11. The model is further rotated 45° counterclockwise about the z-axis (azimuth angle) and 45° counterclockwise about the x-axis (dip angle).
We use the exact phase velocity in orthorhombic media (Tsvankin, 1997):
V2(x,k)|k|2=2a2/3b3cos(υ3+n2π3)a3,
(35)
where n=0 corresponds to the P-wave and n=1, 2 corresponds to S-waves, and
υ=arccos(2(a/3)3ab/3+c2((a2/3b)/3)3)(0υπ),
(36)
a=(G11+G22+G33),
(37)
b=G11G22+G11G33+G22G33G122G132G232,
(38)
c=G11G232+G22G132+G33G122G11G22G332G12G13G23,
(39)
G11=c11k^x2+c66k^y2+c55k^z2,
(40)
G22=c66k^x2+c22k^y2+c44k^z2,
(41)
G33=c55k^x2+c44k^y2+c33k^z2,
(42)
G12=(c12+c66)k^xk^y,
(43)
G13=(c13+c55)k^xk^z,
(44)
and
G23=(c23+c44)k^yk^z.
(45)
To incorporate tilting into the orthorhombic anisotropy, we replace the original wavenumber components kx, ky, and kz with k^x, k^y, and k^z, which are wavenumbers evaluated in the rotated coordinate system aligned with the symmetry axis:
k^x=kxcosϕ+kysinϕ,k^y=kxsinϕcosθ+kycosϕcosθ+kzsinθ,k^z=kxsinϕsinθkycosϕsinθ+kzcosθ,
(46)
where ϕ is the azimuth angle representing horizontal rotation (the angle between the original x-axis and the rotated one), and θ is the dip angle measured from vertical. Figure 15 demonstrates the wavefield snapshots taken at t=1s for three wave modes: the quasi-P-wave and the coupled quasi S-waves. Note that the quasi S-waves are propagated separately using solutions from equation 35 and then summed together because the two modes do not decouple easily in an orthorhombic medium.

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

In this appendix, we prove the unconditional stability of the one-step wave extrapolation linear operator u=Lf in 1D isotropic media defined by
u(x)=ξze2πi(xξ+V(x)|ξ|t)f^(ξ)forf(x)l2[0,1],
(A-1)
where 2πξ=k, f(x) is assumed to have a periodic boundary condition, and f^(ξ) is the Fourier transform of f(x) as defined by equation 3. We treat ξ as discrete and x as continuous for ease of derivation. Our argument can be viewed as a discrete version of the standard stationary phase method in the study of pseudodifferential operators (Stein, 1993; Grigis and Sjöstrand, 1994). To show that operator L is stable, a sufficient condition is that L221+Ct, where C is a bounded constant. From equation A-1, we observe that operator L is the composition of two operators L=AF, where F is the inverse Fourier transform and A is the operator defined by
(Aω)(z)=ξze2πi(xξ+V(x)|ξ|t)ω(ξ)forωl2(z).
(A-2)
Let us consider AA:l2l2, where A corresponds to a matrix with (x,ξ) entry given by A(x,ξ)=e2πi(xξ+V(x)|ξ|t) and A corresponds to a matrix with (η,x) entry given by A(η,x)=e2πi(xη+V(x)|η|t). The AA represents a matrix with (η,ξ) entry given by
AA(η,ξ)=e2πi[x(ξη)+V(x)(|ξ||η|)t]dx.
(A-3)
To bound the l2l2-norm of AA, we estimate the (η,ξ) entry of AA. For ξ=η, we have AA(η,ξ)=1. For ξη
x(ξη)+V(x)(|ξ||η|)t=x(ξη)+V(x)α(ξη)t,
(A-4)
with α=(|ξ||η|)/(ξη). Clearly, |α|1. Then, AA can be expressed as
AA(η,ξ)=e2πi(ξη)[x+V(x)αt]dx.
(A-5)
For sufficiently small t, x+V(x)αt satisfies
12x[x+V(x)αt]32.
(A-6)
Equation A-5 can be expressed as
AA(η,ξ)=e2πi(ξη)[x+V(x)αt]1+V(x)αt1+V(x)αtdx,=11+V(x)αte2πi(ξη)[x+V(x)αt]d[x+V(x)αt].
(A-7)
Let us define y=x+V(x)αt. From equation A-6, it is clear that the map xy is one to one. Substituting y into equation A-7 gives
AA(η,ξ)=11+V(x(y))αte2πi(ξη)ydy,
(A-8)
which is the inverse Fourier transform of 11+V(x(y))αt. When t is small,
AA(η,ξ)=[1V(x(y)αt+O(t2))]e2πi(ξη)ydy,=αtV(x(y))e2πi(ξη)ydy+O(t2).
(A-9)
To evaluate the norm of the integration term in the last equation, we perform integration by part for k-times and apply periodic boundary condition:
|V(x(y))e2πi(ξη)ydy|1(2π|ξη|)ke2πi(ξη)y|ykV(x(y))|dy.
(A-10)
Assuming sufficient smoothness on V(x), we have for ξη
|(AAI)(η,ξ)|Ct|ξη|k,
(A-11)
for a constant C. To estimate the l2l2-norm of AA, we make use of the following lemma, which can be derived from direct calculation:

We now apply the above lemma to AAI. When kd, where d is the spatial dimension, we have ξηC|ξη|k bounded. Therefore, for sufficiently smooth V, we have AAI22Ct, for sufficiently small t. Hence, AA221+Ct and A221+Ct. Because L=AF and F as the Fourier transform is an isometry, we have L221+Ct.

When performing wave extrapolation, fix a final time T and propagate T/dt steps, the operator is stable because
LT/t22(1+Ct)T/teCT.
(A-12)
Freely available online through the SEG open-access option.