In the common-shot gathers, the reflection and the diffraction generated by a point source share similar kinematic characteristics, which are difficult to be distinguished (Landa et al., 1987). While when the incident wave is a plane wave, the travel time curve of the reflected wave can be written as: t R = cos ⁡ θ ⁡ sin ( θ − 2 ϕ ) v ⁡ cos ( θ − ϕ ) ⋅ x + 2 h v [ cos ⁡ θ ⁡ cos ( θ − ϕ ) − c o s ϕ ⁡ cos ( θ − ϕ ) − cos ⁡ θ cos ( θ − ϕ ) ] (1) where t R denotes the travel time and v denotes the medium velocity. The specific meaning of θ and ϕ is shown in Figure 1.

The travel time curve of the diffracted wave can be expressed as: ( t D − h / v ) 2 h 2 ⁡ cos 2 ⁡ θ / v 2 − ( x − h ⁡ sin ⁡ θ ) 2 h 2 ⁡ cos 2 ⁡ θ = 1 (2)

According to Eqs. 1, 2 and Figure 1, the seismic response of the reflection is linear while the response of the diffraction is hyperbolic. Therefore, we can separate the reflection and the diffraction based on the difference of kinematic properties in plane-wave gathers. Firstly, we convert the common-shot gathers into common ray-parameter gathers using tau-p transform (Kappus et al., 1990; Zhang et al., 2005). The tau-p transform equation can be expressed as: P ( τ , p ) = ∫ d ( t , x ) d x = ∫ d ( τ + p x , x ) d x (3) where d ( t , x ) denotes the seismic shot records in the time-offset t − x domain, and P ( τ , p ) is the so-called plane-wave gathers with ray parameter p and time axis τ . Then we use the PWD wave filter proposed by Fomel (2002) to estimate the local dip angle of the diffracted wave and separate it in the tau-p domain. After that, we use the inverse tau-p transform to get the diffraction wavefield in the time-offset domain. The inverse tau-p transform equation is expressed as: d D ( t , x ) = ∫ P ( τ , p ) d p = ∫ P ( τ = t − p x , p ) d p (4) where d D ( t , x ) denotes the shot records of the diffracted wave.

In Tarim Basin, where the fractured-vuggy carbonate reservoirs are extensively developed, the surface is often covered by desert and gravel layers with low velocity. In forward modeling and migration, a fine spatial grid must be used to avoid numerical dispersion. However, for the deep region with high velocity, using fine grid is a waste of computing resources. We use an adaptive sampling strategy to solve the above problem.

As shown in Figure 2, the horizontal axis denotes depth and the vertical axis represents vertical grid interval. The black dots on the x-axis denote initial vertical grid points. The black solid line denotes initial vertical grid interval, which is fixed. According to the medium velocity and the dominant frequency of the source, we use Eq. 5 to recalculate the optimal vertical grid interval (The red solid line in Figure 2): d z ( z ) = v min ( z ) 10 f m (5) where v min ( z ) denotes the minimum velocity along the depth axis z, f m is the dominant frequency of the source wavelet, and d z ( z ) denotes the optimal vertical grid spacing.

We use a rectangular sampling method (Wang et al., 2020) to resample the initial migration model. Firstly, we set a small trial step from z = 0 and increase it continuously to get the first grid point η 1 , where η 1 = d z ( η 1 ) . Then we repeat the previous step to get the second grid point η 2 , where η 2 − η 1 = d z ( η 2 ) − d z ( η 1 ) . Finally, we repeat the above process to the max depth and obtain a new model. The vertical grid interval of the new model is irregular and it varies with velocity. Note that there is a mapping relationship between the initial model and the new model, we use two different coordinate systems to express it: x = x ( ξ , η ) (6) z = z ( ξ , η ) (7) where x and z are the coordinate variables in coordinate system A ( x , z ) , ξ and η are the coordinate variables in coordinate system B ( ξ , η ) . The initial model is located in A ( x , z ) and the new model is located in B ( ξ , η ) . The derivation process of the mapping relationship is shown in Supplementary Appendix SA.

The acoustic wave equation in 2D heterogeneous isotropic medium in coordinate system A ( x , z ) can be written as: 1 v 2 ∂ 2 u ∂ t 2 = ρ ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ z 2 ) + f s ( t ) (8) where u denotes acoustic pressure field, t denotes time, ρ denotes medium density, v denotes acoustic velocity, and f s ( t ) denotes source function.

Substitute Eq. A-9 and A-10 into Eq. 8, we get the expression of Eq. 8 in coordinate system B ( ξ , η ) : 1 v 2 ∂ 2 u ∂ t 2 = ρ ( ∂ 2 u ∂ ξ 2 + ∂ 2 u ∂ η 2 1 z η 2 − ∂ u ∂ η z η η z η 3 ) + f s ( t ) (9) where z η = ∂ z / ∂ η and z η η = ∂ 2 z / ∂ η 2 they can be solved by finite-difference method.

Compared to Eq. 8, one term is added to Eq. 9, which increases the computational cost for one grid point. In fact, the omission of the extra term has little effect on the final imaging results. Therefore, Eq. 9 is further simplified as: 1 v 2 ∂ 2 u ∂ t 2 = ρ ( ∂ 2 u ∂ ξ 2 + ∂ 2 u ∂ η 2 1 z η 2 ) + f s ( t ) (10)

LSRTM is considered to be a true amplitude imaging method, which can improve imaging resolution, reduce artifacts, and meet the demand for lithological interpretation.

As shown in Eq. 11, the seismic wavefield satisfies the superposition principle: u = u 0 + u s (11) where u denotes total wavefield, u 0 denotes background wavefield, and u s denotes perturbation wavefield.

Similarly, the velocity model can be regarded as the superposition of a perturbation model and a smoothed background model, 1 v 2 = 1 v 0 2 + 1 Δ v 2 (12)

u 0 obeys Eq. 10: 1 v 2 ∂ 2 u 0 ∂ t 2 = ρ ( ∂ 2 u 0 ∂ ξ 2 + ∂ 2 u 0 ∂ η 2 1 z η 2 ) + f s ( t ) (13)

Substitute Eq. 11 and Eq. 12 into Eq. 10, subtract Eq. 13, and apply Born approximation (Dai et al., 2012), we can obtain the control equation of 1 v 0 2 ∂ 2 u s ∂ t 2 = ρ ( ∂ 2 u s ∂ ξ 2 + ∂ 2 u s ∂ η 2 1 z η 2 ) − 1 Δ v 2 ∂ 2 u 0 ∂ t 2 (14)

Eq. 14 is the Born (linearized) forwarding modeling equation, and it can be rewritten as a matrix: d s = L m (15) where L is the Born forwarding modeling operator, m denotes model parameter, and d s denotes the Born-modeled data. We can get the migration image by using the the RTM operator: m m i g = L T d s (16) where m m i g represents the imaging result. L T denotes the RTM operator, which is the conjugate transpose of L .

The goal of LSM theory is to minimize the objective function, which is defined as: J ( m ) = 1 2 ‖ L m − d o b s ‖ 2 2 (17) where J ( m ) is the objective function, d o b s denotes the observed data, and ‖ ⋅ ‖ 2 2 is the L2-norm of a vector.

We use a conjugate-gradient algorithm (Dai et al., 2012) to solve the model parameter m , which is formulated as: g k + 1 = L ∗ [ L m k − d o b s ] β k + 1 = [ g k + 1 ] ∗ [ C g k + 1 ] [ g k ] [ C g k ] z k + 1 = C g k + 1 + β k + 1 z k a k + 1 = [ z k + 1 ] ∗ [ g k + 1 ] [ L z k + 1 ] ∗ [ L z k + 1 ] m k + 1 = m k − a k + 1 z k + 1 (18) where k and * denote iteration index and conjugate transpose of a matrix, g k and β k represent gradient and correction coefficient, z k and α k denote conjugate gradient direction and step length, and C denotes precondition operator.

In summary, the implementation of the AS-D-LSRTM method contains four main steps. The first step is to separate diffraction. The second step is to resample the initial model. The third step is to obtain the image after several iterations by LSRTM. Finally, we use linear interpolation to convert the final image from the coordinate system B ( ξ , η ) to coordinate system A ( x , z ) . Figure 3 shows the complete workflow of AS-D-LSRTM.

The initial fault model is shown in Figure 4A, which has a low-velocity surface layer with velocity of 2000 m/s. It is discretized on a 801 × 301 grid a with grid spacing d x = d z = 8 m. The time sampling interval is 0.5 ms and the recording time is 1.5 s. A total of 51 sources are distributed laterally from 2.4 to 4 km. The shot interval is 32 m. There are 601 receivers of each shot and the receiver interval is 8 m. We use a Ricker wavelet as the source, and its dominant frequency is 25 Hz. To improve the computational efficiency, we resample the initial model using the adaptive sampling method. Figure 4B shows the resampled model with 151 vertical grid points. From Figure 4B we can see that the deep layers with high velocity are compressed. Figure 5 shows the grid interval comparison of the two models. The black line in Figure 5 denotes initial vertical grid interval, the blue line denotes theoretical value calculated by Eq. 5, and the red line shows the vertical grid interval of the resampled model.

We use the tau-p transform to produce 501 plane-wave time sections with ray parameters ranging from −0.5 to 0.5 ms/m. Figure 6A displays a section with ray parameter p = 0  ms/m. In Figure 6A, the reflection is linear (see the black arrows) and the diffraction is hyperbolic (see the red arrows). Figure 6B shows the estimated dip field of Figure 6A. Based on the dip angle field shown in Figure 6B, we use the PWD filter method to separate the diffraction. Figure 6C shows the separated diffraction with ray parameter p = 0  ms/m. We can see from Figure 6C that the reflection is removed while the diffraction is preserved. Figure 7A shows the point source response of the full wavefield. In Figure 7A, the diffraction energy generated by the inclined faults (the red arrow) is much weaker than that of the layer interface (the black arrow). Figure 7B shows the response of the diffraction wavefield. In Figure 7B, the reflection is suppressed and the diffraction (the blue arrow) is preserved.

Figures 8A, 8B show the images of full wavefield LSRTM (F-LSRTM) and D-LSRTM after 30 iterations, respectively. In Figure 8B, we can see that the acquisition footprints (the red arrow) and the reflection generated by the surface layer (the blue arrow) are almost eliminated. Moreover, the resolution of the fault in Figure 8B is higher than that in Figure 8A (see the black arrows). Figure 8C shows the AS-D-LSRTM image after 30 iterations, and Figure 8D shows the final result after linear interpolation. In Figure 8B, the artifacts and the diffraction are enhanced simultaneously when the reflection is removed. Nevertheless, the diffraction-based LSRTM results are still focused. Figure 9 displays the normalized residual convergence curves of F-LSRTM, D-LSRTM, and AS-D-LSRTM. The curves show that the convergence rate of D-LSRTM and AS-D-LSRTM is faster than that of F-LSRTM. Moreover, the normalized residual of D-LSRTM and AS-D-LSRTM decreases quickly before the 5th iterations and gradually converges to a same value. Figure 10 shows the normalized computing time and memory footprint of these methods. Compared with D-LSRTM, AS-D-LSRTM can save 44% of computing cost and half of memory requirement without losing accuracy. From the above analysis, we conclude that AS-D-LSRTM is superior to F-LSRTM and D-LSRTM.

We use a layered model to further verify the effectiveness of the proposed method. Figure 11A shows the initial model, which is discretized on a 801 × 401 grid with a grid interval of d x = d z = 5 m. The model has two layers with velocity of 2.5 km/s and 4 km/s, respectively. There are three small-scale vugs and four fractures distributed at the layer interface with velocity of 3 km/s. The vugs are square with side lengths of 20, 40, and 60 m, respectively. The width of the fractures is 20 m and the lengths of them are 100, 200, 300, and 400 m, respectively. Figure 11B shows the resampled model ( f m in Eq. 5 is 30), which has 221 vertical grid points in total. In total, 41 sources are distributed laterally from 1.5 to 2.5 km. Each shot has 601 receivers. The shot interval is 25 m and the receiver interval is 5 m. The time step is 0.3 ms and the recording time is 1.5 s.

In this tests, we apply the adaptive sampling method to F-LSRTM (denoted by AS-F-LSRTM) and D-LSRTM. Figure 12 shows the images after 30 iterations using different LSRTM methods. In Figures 12A,C, the vugs and fractures are imaged successfully, and the bottom of the fractures is shown as a string of beads (Lu et al., 2015). In Figures 12B,D, the reflection are restrained and the energy of the diffraction is enhanced.

Figure 13 displays the normalized residual convergence curves of the four methods. The result shows that the convergence rate of D-LSRTM and AS-D-LSRTM is faster than that of F-LSRTM and F-AS-LSRTM. We can see that the convergence curves of D-LSRTM and AS-D-LSRTM are almost coincident after the 16th iteration. Figure 14 shows the normalized computing cost and storage requirement of the four methods. From Figures 12–14, we conclude that the AS-D-LSRTM algorithm is helpful to identify small-scale caves and fractures, and it can greatly save computing resources.

The fault-karst carbonate reservoirs are usually stripped distribution along the strike-slip fault zone, especially in South Tahe area of Tarim Basin (Ding et al., 2020). According to stress state, the strike-slip fault can be divided into three structural styles, which are translation, extrusion and extension. The oil and gas resources are more abundant in the extension and extrusion section (Cheng et al., 2020). Figure 15 shows the distribution characteristics of the fault-karst reservoirs. In the extension section (left), large-scale caves are developed near the trunk fracture zone and small-scale vugs are distributed along the branch fracture. In the extrusion section (right), the caves and vugs are mainly distributed along the trunk fracture. According to the characteristics of the fault-karst carbonate reservoirs formed by strike-slip extension, we build an ultra-deep fault-karst model, as shown in Figure 16A. The model includes 1101 × 1101 grid points with grid spacing of 8 m. The maximum depth of the fault-karst structure is more than 8 km. Figure 16B shows the new model after adaptive sampling. The number of vertical grid points reduces from 1101 to 510.

In numerical tests, a total of 61 sources are distributed laterally from 3.2 to 5.6 km. Each shot has 801 receivers. The shot interval is 40 m and the receiver interval is 8 m. The time step is 0.5 ms and the record time is 4 s. The source is a Ricker wavelet and its dominant frequency is 25 Hz. Figures 17A,B show AS-F-LSRTM and AS-D-LSRTM images after 20 iterations, respectively. In Figure 17A, some small-scale caves and vugs are shown as strong reflection spots (see the red arrows). However, the reflection energy generated by large-scale caves near the trunk fracture is much weaker (see the black arrows). In addition, the reflection-based LSRTM images cannot correctly describe the distribution and the structural characteristics of the fault-karst reservoirs. For example, the reflection energy of the top of the reservoirs is extremely weak (see the blue arrows in Figure 17A). This phenomenon can mislead the prospectors and lead them to miss productive reservoirs. In Figure 17B, the diffraction of the trunk fracture is much stronger than that in Figure 17A (see the black arrows). Moreover, the top of the reservoirs is easier to identify in Figure 17B (see the blue arrows). We conclude that AS-D-LSRTM can enhance the diffraction and improve the imaging resolution of the fault-karst carbonate reservoirs.

The LSRTM algorithms relies heavily on accurate migration velocity model and high-quality data. In this paper, the proposed AS-D-LSRTM method is also tested on a 2D field dataset to further verify its effectiveness. Figure 18A shows the initial migration velocity model, which includes 2001 × 401 grid points with a grid spacing of 10 m. Figure 18B shows the variable-grid model after adaptive sampling (n in Eq. 5 is 10), which is discretized on a 401 × 284 grid. We use the Ricker wavelet with dominant frequency of 20 Hz as the source. Total shots are 100 with irregular distribution and the record time is 4 s. Each shot has 204 receivers and the receiver interval is 50 m. The time sampling interval is 0.5 ms.

Figure 19A shows a single shot record. We apply the PWD filter method to the full-wavefield data and obtain its diffraction component, as shown in Figure 19B. Figures 20A,B show the AS-F-LSRTM and AS-D-LSRTM images after 5 iterations, respectively. In Figure 20A, the diffraction generated by the fault (see the red arrow) is muck weaker than the reflection produced by the horizontal layers (see the black arrow). From the image result of AS-D-LSRTM (Figure 20B), we can see that the reflection energy (see the red arrow) is almost eliminated while the diffraction is enhanced (see the black arrow). In addition, about 74% of computing time and 71% of memory space are saved after using the adaptive sampling method. We conclude that our method is still effective for field data.