Reverse time migration (RTM) was proposed in the 1980s (Baysal, 1983; McMechan, 1983; Whitmore, 1983), which is based on a two-way wave equation and applies zero-lag cross-correlation imaging conditions to realize imaging. Theoretically, RTM can adapt to any complex velocity model without dip limitations and image nearly all kinds of waves, including refractions, prismatic waves, diffractions, and multiples, so it is considered to be the most accurate imaging algorithm and has been widely used in the field data processing (Sun et al., 2015; Oh et al., 2018; Qu et al., 2020; Fee et al., 2021). However, due to using the two-way wave equation to implement wavefield continuation, the backward reflection will occur when the seismic wave propagates to the reflection interface. The conventional zero-lag cross-correlation imaging conditions directly use all forward- and back-propagated wavefields to form subsurface images (Du et al., 2013; Chen and He, 2014; Fei et al., 2015), which could inevitably produce a lot of low-wavenumber migration artifacts.

At present, three main methods can be used to reduce the migration artifacts in RTM. The first one is a backward reflection suppression method, which usually employs the nonreflecting acoustic equation to imaging. Baysal (1984) has first proposed a nonreflecting acoustic equation based on an assumption of constant wave impedance, which can significantly suppress the backward reflection of the vertical incident seismic waves. Song (2005) has improved the nonreflecting acoustic equation to enhance the suppression effect of backward reflection. However, in general, the backward reflection suppression effect is not ideal, and it is difficult to achieve the purpose of effectively eliminating migration artifacts. The second one is the filtering method. Mulder and Plessix (2004) have directly used high-pass filtering to denoise the imaging section. Zhang and Sun (2009) have applied Laplacian filtering to filter the results of RTM. However, this kind of method has problems such as difficulty in determining the threshold, damage to the effective signal, and incomplete noises removal. The third one is to modify the imaging conditions. There are two kinds of methods used to modify the imaging conditions usually. One is the angle weighting method proposed by Yoon and Marfurt (2006), which can effectively remove the large-angle migration artifacts by introducing an attenuation factor related to the reflection angle into the imaging conditions. The other is the wavefield decomposition method, which decomposes the source and receiver wavefields into going wavefields in different directions and then extracts the effective wavefield components to form images to achieve the accurate imaging of underground structures. Some scholars (Liu et al., 2011; Fei et al., 2015; Wang et al., 2016) have successively applied the Hilbert transform to realize wavefield decomposition and obtained high-precision imaging sections. However, when the Hilbert transform is applied to wavefield decomposition, a certain amount of the computational cost is required. Chen and He (2014) have used the Poynting vector to decompose the source and receiver wavefields in the four directions of up, down, left, and right, which can greatly improve the suppression effect of low-wavenumber noises with small additional computation cost. Therefore, this method has been widely used (Wang and He, 2017; Liu, 2019; Li and He, 2020; Wang et al., 2021). However, there are also two problems in wavefield decomposition based on the Poynting vector method. First, it is not accurate enough for the Poynting vector method to indicate all directions of seismic wave propagation (Du et al., 2012; Zhang, 2014; Duan and Sava, 2015; Li and He, 2020) and the second is that there are always some singularities in the Poynting vector.

The optical flow method was first proposed to solve the motion information problem of objects between adjacent frames (Horn and Schunck, 1981; Lucas and Kanade, 1981), and then it was introduced to RTM (Hu et al., 2014; Zhang, 2014; Gong et al., 2016; Zhang et al., 2018; Wu et al., 2021). Compared with the Poynting vector, the optical flow vector is a more accurate vector that is closer to the real wavefield propagation direction. Moreover, there is no singularity in the optical flow vector. In this article, the optical flow vector method is introduced into RTM to decompose wavefields and calculate the reflection angle of each imaging point underground. Based on the optical flow vector method, both source and receiver wavefields can be decomposed accurately and the accurate reflection angle of each imaging point underground can be obtained; then, by the introduction of an attenuation factor related to the reflection angles, the angle-weighted RTM with wavefield decomposition based on the optical flow vector is implemented, which greatly improves RTM imaging.

In the next section, we review the wavefield continuation of RTM based on the acoustic wave equation. In Wavefield Decomposition Based on the Optical Flow Vector, the wavefield decomposition based on the optical flow vector method is introduced and some tests are given to compare the effects of wavefield decomposition for the Poynting vector method and the optical flow vector method. In Angle-Weighted RTM Imaging Based on the Optical Flow Vector, we show how to calculate the reflection angle of each imaging point underground based on the optical flow vector method and how to produce the final RTM image using an attenuation factor related to the reflection angles. In Numerical Tests on the Marmousi Model and Field Marine Seismic Data Imaging, we present some tests to show the imaging effect of the method developed in the article. We end with some concluding remarks in Conclusion.

The first-order stress-velocity acoustic wave equation in a two-dimensional isotropic medium can be expressed as follows: { ∂ v x ∂ t = 1 ρ ∂ p ∂ x ∂ v z ∂ t = 1 ρ ∂ p ∂ z ∂ p ∂ t = ρ v 2 ( ∂ v x ∂ x + ∂ v z ∂ z ) , (1) where x and z represent the space coordinates, respectively; v _x and v _z denote the practical vibration velocity in the x and z direction, respectively; t is the time; ρ signifies the density; v represents the velocity of the acoustic wave; p denotes the stress.

We use staggered grids to discretize Eq. 1 by finite-difference for realizing forward wavefield continuation and reverse time wavefield continuation. Taking forward continuation as an example, the high-order difference schemes of Eq. 1 can be written as follows: { v x k ( i + 1 / 2 , j ) = v x k − 1 ( i + 1 / 2 , j ) + Δ t ρ ( i , j ) Δ x ∑ m = 1 N C m [ p k − 1 / 2 ( i + m , j ) − p k − 1 / 2 ( i − m + 1 , j ) ] v z k ( i , j + 1 / 2 ) = v z k − 1 ( i , j + 1 / 2 ) + Δ t ρ ( i , j ) Δ z ∑ m = 1 N C m [ p k − 1 / 2 ( i , j + m ) − p k − 1 / 2 ( i , j − m + 1 ) ] p k + 1 / 2 ( i , j ) = p k − 1 / 2 ( i , j ) + ρ ( i , j ) v 2 ( i , j ) Δ t Δ x ∑ m = 1 N C m [ v x k ( i + m − 1 / 2 , j ) − v x k ( i − m + 1 / 2 , j ) ] + ρ ( i , j ) v 2 ( i , j ) Δ t Δ z ∑ m = 1 N C m [ v z k ( i , j + m − 1 / 2 ) − v z k ( i , j − m + 1 / 2 ) ] , where k represents the temporal discrete point number, i and j denote the spatial discrete point numbers in the x and z direction, respectively. Δt is the time discrete step; Δx and Δz are the spatial discrete steps in the x and z directions, respectively. N denotes half of the accuracy of spatial difference, and C _m is the difference coefficients.

In wavefield continuation based on the finite-difference method, artificial boundaries have been used in practice to suppress boundary reflection. To eliminate the boundary reflection, the perfectly matched layer (PML) method is used here. PML boundary algorithm has been widely studied (Berenger, 1994; Collino and Tsogka, 2001; Zhang and Shen, 2010), so we do not discuss it in detail.

The Poynting vector, also known as the energy flux density vector, was first applied in the field of electromagnetic computing (Poynting, 1884). Now, it has become a common algorithm used to indicate the propagation direction of wavefields in seismic wavefield calculation (Yoon and Marfurt, 2006; Tang et al., 2017).

The Poynting vector of the first-order stress-velocity acoustic wave equation can be expressed as follows: P y  = - ∇ u ∂ u ∂ t = ( − p v x , − p v z ) = ( P x y , P z y ) , (3) where u represents the wavefields, P^y _z and P ^y _z are the horizontal and vertical components of the Poynting vector, respectively, and ∇ denotes the gradient operator. We can obtain, using the Poynting vector, the up-, down-, left-, and right-going wavefields after wavefield decomposition. Taking the wavefield decomposition of source wavefields as an example, the specific expression can be represented as follows: S u ( x , z , t ) = { S ( x , z , t )     i f   P z y < 0    0 S d ( x , z , t ) = { S ( x , z , t )     i f   P z y ≥ 0    0 S l ( x , z , t ) = { S ( x , z , t )     i f   P x y < 0   0 S r ( x , z , t ) = { S ( x , z , t )     i f   P x y ≥ 0     0 , (4) where S _u (x, z, t), S _d (x, z, t), S _l (x, z, t), and S _r (x, z, t) are the up-, down-, left- and right-going source wavefields, respectively. It can be seen from Eq. 3 that the calculation of the Poynting vector is composed of the product of the time derivative and the space derivative of the wavefield. When the time derivative or the space derivative is zero, the Poynting vector must be zero too, which causes instability. Furthermore, Zhang et al. (2014) have pointed out that the Poynting vector itself is difficult to indicate the propagation direction of the wavefield with high accuracy.

The optical flow vector is a vector that is obtained by several iterations and can indicate the propagation direction of the wavefield stably and accurately. Therefore, we introduce the optical flow vector into the wavefield decomposition process of RTM. In the two-dimensional RTM, the fundamental assumption for the optical flow problem is that the wavefield u at a spatial point (x, z) is continuous for very small variations in space (dx and dz) and time (dt), and its expression is as follows: u ( x + d x , z + d z , t + d t ) = u ( x , z , t ) , (5) where u denotes the wavefields, x and z represent the space coordinates, respectively, and t is time. We use the Taylor formula to expand u (x + dx, z + dz, t + dt) and discard higher-order terms above the second order and obtain u x P x o + u z P z o + u t = 0 , (6) where u _x and u _z are the spatial derivatives of wavefields, u _t is the time derivatives of wavefields, and P^o _x and P ^o _z are the orthogonal (x and z) components of the optical flow vector, respectively. With two unknowns (P^o _x and P ^o _z ) and only one Eq. 6, the problem is ill-posed and the solution is nonunique. To address this underdetermined problem (Eq. 6), the regularization terms of global smooth constraints are introduced by requiring that neighboring points have similar flow directions as that at a central target point. Therefore, we construct the following misfit function: E = ∬ [ ( u x P x o + u z P z o + u t ) 2 + α 2 C ] d x d z , (7) where C is the regularization terms, which can be written as follows: C = ∇ 2 P x o + ∇ 2 P z o , (8) where ∇² is the Laplacian operator and α is a weighting factor of the regularization term, generally taken as 1. Equation 7 can be solved using an iterative least-squares approach, in which the update parameters are computed as follows: { ( P x o ) n + 1 = ( P ¯ x o ) n − u x [ u x ( P ¯ x o ) n + u z ( P ¯ z o ) n + u t ] α 2 + u x 2 + u z 2 ( P z o ) n + 1 = ( P ¯ z o ) n − u z [ u x ( P ¯ x o ) n + u z ( P ¯ z o ) n + u t ] α 2 + u x 2 + u z 2 , (9) where ‾P^o _x and‾ P ^o _z are the local average values of the horizontal and vertical components of the optical flow vector, respectively, and n is the number of iterations. From Eqs 7, 8, and 9, it can be seen that the optical flow vector is no longer a simple vector generated by multiplying the time derivative and the spatial derivative of the wavefield directly, but an accurate vector obtained by several iterative operations based on an initial optical flow vector, so it is closer to the real propagation direction information of the wavefield. Besides, because of the addition of the regularization term in the calculation process of the optical flow vector, only when the time derivative and space derivative are both zero, the optical flow vector is zero, which avoids the instability effectively in the calculation process of the optical flow vector.

The feasibility and accuracy of the method are first evaluated using a two-layer velocity model (as shown in Figure 1). The size of the homogeneous medium model is 1,500 m in length and 1,500 m in depth. The velocity of the first layer is 2,500 m/s and the second layer is 3,000 m/s. A Ricker wavelet with a dominant frequency of 30 Hz is used as the source, which is excited at (750 m, 0 m). The grid interval in the x and z directions is 5 m. The finite-difference accuracy of wavefield continuation is tenth order in space. The time sampling step is 0.5 ms and the maximum recording time is 0.8 s. Figure 2 illustrates the source wavefield snapshot at 0.3 s. Figures 3A,B, respectively, show the wavefield direction near the reflection interface (indicated by the red box in Figure 2) calculated using the Poynting vector and the optical flow vector. Figures 4A,B contain the horizontal components of the Poynting vector and the optical flow vector at this time, respectively. The left-going wavefield obtained by wavefield decomposition based on the Poynting vector and the optical flow vector are plotted in Figures 5A,B.

From Figures 3A,B (indicated by the red circle) and Figures 4A,B (indicated by the red arrow), it can be seen that accurately indicating the propagation direction of the wavefield using the Poynting vector is challenging and singular values are prone to appear, whereas the optical flow vector is smoother and the instability phenomenon is avoided effectively. Comparing Figures 5A,B, we can see that for the wavefield decomposition achieved based on the Poynting vector, some other wavefield components as indicated by the arrow appear because the Poynting vector calculation is inaccurate and unstable, whereas the optical flow vector does not generate other wavefield components, so the decomposed wavefield is more accurate.

Based on the optical flow vector, we use Eq. 10 to decompose the source and receiver wavefields in the four directions of up, down, left, and right: S u ( x , z , t ) = { S ( x , z , t )     i f   P z o < 0    0 S d ( x , z , t ) = { S ( x , z , t )     i f   P z o ≥ 0    0 S l ( x , z , t ) = { S ( x , z , t )     i f   P x o < 0    0 S r ( x , z , t ) = { S ( x , z , t )     i f   P x o ≥ 0    0 R u ( x , z , t ) = { R ( x , z , t )     i f   P z o < 0    0 R d ( x , z , t ) = { R ( x , z , t )     i f   P z o ≥ 0    0 R l ( x , z , t ) = { R ( x , z , t )     i f   P x o < 0    0 R r ( x , z , t ) = { R ( x , z , t )     i f   P x o ≥ 0    0 , (10) where S _u (x, z, t), S _d (x, z, t), S _l (x, z, t), and S _r (x, z, t) denote the up-, down-, left-, and right-going source wavefields, respectively, and R _u (x, z, t), R _d (x, z, t), R _l (x, z, t), and R _r (x, z, t) represent the up-, down-, left-, and right-going wavefields of receivers, respectively.

The decomposed wavefields of both sources and receivers in the opposite direction are selected for imaging separately to avoid migration artifacts (Chen and He, 2014), using the following: I o p ( x , z ) = ∑ t S u ( x , z , t ) R d ( x , z , t ) + ∑ t S d ( x , z , t ) R u ( x , z , t ) ∑ t S 2 ( x , z , t ) + + ∑ t S l ( x , z , t ) R r ( x , z , t ) + ∑ t S r ( x , z , t ) R l ( x , z , t ) ∑ t S 2 ( x , z , t ) , (11) where I _op (x, z) is the optimal stack section. However, in fact, the migration artifacts are mainly distributed in the large-angle region (Yoon and Marfurt, 2006; Zhang et al., 2014; Zhang et al., 2020). Although the method of selecting the wavefields in the opposite direction for imaging can reduce the migration artifacts of 180° or close to 180°, the suppression effect on the migration artifacts coming from other large-angle regions is weak. Moreover, there is also a risk of losing effective information if only the wavefields in the opposite direction are selected for imaging. To address this problem, the reflection angle of each imaging point underground is first calculated based on the optical flow vector. The calculation formula can be defined as follows: θ = π − a r cos ( P S o ( x , z , t ) · P R o ( x , z , t ) | P S o ( x , z , t ) | ∗ | P R o ( x , z , t ) | ) , (12) where θ is the reflection angle of each imaging point and P^o _S and P^o _R are the optical flow vectors of the source and receiver wavefields, respectively. Then, an attenuation factor related to the reflection angle is introduced as a weight to generate the optimal stack section, and the final imaging result is obtained according to Eq. 13: I ( x , z ) = ( ∑ t S u ( x , z , t ) R d ( x , z , t ) + ∑ t S d ( x , z , t ) R u ( x , z , t ) ∑ t S 2 ( x , z , t ) + + ∑ t S l ( x , z , t ) R r ( x , z , t ) + ∑ t S r ( x , z , t ) R l ( x , z , t ) ∑ t S 2 ( x , z , t ) ) w ( θ ) , (13) where I (x, z) is the final imaging result of angle-weighted RTM with wavefield decomposition based on the optical flow vector and w(θ) is the attenuation function, and we choose a cosine-type function as the attenuation function.

The two-layer velocity model in Wavefield Decomposition Based on the Optical Flow Vector is used to test the effect of the angle-weighted imaging method. Figure 6A shows the result of conventional RTM, Figure 6B illustrates the result of RTM with wavefield decomposition based on the Poynting vector, Figure 6C contains the result of RTM with wavefield decomposition based on the optical flow vector, and Figure 6D is the result of angle-weighted RTM with wavefield decomposition based on the optical flow vector.

There are obvious migration artifacts in the conventional RTM imaging result in Figure 6A. As shown in Figures 6B–D, we can see that most of the migration artifacts are eliminated in the result of RTM with wavefield decomposition based on the Poynting vector. However, due to inaccurate wavefield decomposition, there are still some noises remaining, and as a result of RTM with wavefield decomposition based on the optical flow vector, the migration artifacts are further suppressed. Moreover, the migration artifacts are basically completely suppressed, and the effective structural imaging is highlighted by performing angle weighting processing on the optimal stack section. Therefore, the angle-weighted RTM with wavefield decomposition based on the optical flow vector can produce the accurate imaging of underground structures.

A region of the Marmousi-II model (as shown in Figure 7) is used to test the imaging accuracy of our method. The size of the model is 6,500 m in length and 3,500 m in depth. The grid spacing is 5 m. There are a total of 101 shots and each shot contains 1,300 receivers. The sampling intervals for the shots and receivers are 65 m and 5 m, respectively. The depths of shots and receivers are both 0 m. A Ricker wavelet with a dominant frequency of 30 Hz is used as the source. The time sampling interval is 0.4 ms and the total recording time is 4 s. The finite-difference accuracy of wavefield continuation is eighth order in space. Figure 8A contains the result of conventional RTM, Figure 8B illustrates the result of RTM with wavefield decomposition based on the Poynting vector, Figure 8C shows the result of RTM with wavefield decomposition based on the optical flow vector, and Figure 8D is the result of angle-weighted RTM with wavefield decomposition based on the optical flow vector.

As shown in Figure 8A, the image suffers from low-wavenumber noises. The migration artifacts seriously affect the imaging quality and the real imaging structure is completely concealed. It can be seen from Figures 8B–D that all three methods can significantly suppress migration image noises. However, as shown by the red circle in Figure 8, there are still lots of image noises in Figure 8B, and although the migration artifacts in Figure 8C are further eliminated, a few noises are still left. The noises suppression effect in Figure 8D is the best, the underground structure is the clearest, and the quality of the migration section is greatly improved. Moreover, our attenuation factor puts more weight on the RTM result in the deep part because the reflections generated in the deep part usually have a smaller reflection angle than those in the shallow part for a fixed offset. Therefore, the deep imaging accuracy is further enhanced using our attenuation factor.

A field marine seismic line in the East China sea is selected for the RTM test. The line involves 1,637 shots, among which shots are arranged on the right side, while receivers are on the left side. A total of 648 receivers are allotted for each shot. The interval between shots is 37.5 m and the interval between receivers is 12.5 m. The depths of shots and receivers are both 12.5 m. The minimum offset is 187.5 m and the maximum recording time of shot gather is 8 s. The finite-difference difference accuracy of wavefield continuation is eighth order in space and second order in time. Meanwhile, the time sampling step is 1 ms. Figure 9 shows the seismic record of the 601st shot. Figure 10 illustrates a source wavelet that is extracted from the original data.

Figure 11 shows the velocity model of field data, which is obtained by full waveform inversion. Figure 12 illustrates the RTM sections for field marine seismic data (the part ranging from 10 to 70 km is displayed). Among them, Figure 12A is the result of conventional RTM; Figure 12B illustrates the result of RTM with wavefield decomposition based on the Poynting vector; Figure 12C is the result of RTM with wavefield decomposition based on the optical flow vector; Figure 12D shows the result of angle-weighted RTM with wavefield decomposition based on the optical flow vector. From Figure 12A, we can see that there are lots of low-wavenumber noises in the shallow part, as shown by the black dotted circle, which seriously reduces the imaging quality. It can be seen from Figure 12B that low-wavenumber noises are reduced a lot. In Figure 12C, there are fewer noises than in Figure 12B, and in Figure 12D, there are no obvious noises. Therefore, we can conclude that the method in this article can more effectively eliminate low-wavenumber noises compared to other methods and it is suitable for RTM of real data.

The decomposition of source and receiver wavefields can be accurately implemented using the optical flow vector. Then, the cross-correlation imaging sections of one-way propagation components of the forward- and back-propagated wavefields are optimized and stacked. Furthermore, the reflection angle of each imaging point is calculated based on the optical flow vector, and an attenuation factor related to the reflection angle is used as the weight to give the optimal stack images. The numerical experimental results demonstrate the following:

1) The optical flow vector can be used to decompose the wavefield accurately and stably, and RTM with wavefield decomposition based on the optical flow vector can alleviate the effect of low-wavenumber noises effectively.

The proposed method can be applied to elastic-wave RTM and can be further extended to least-squares RTM.