ORIGINAL RESEARCH article

Front. Earth Sci., 12 May 2025

Sec. Solid Earth Geophysics

Volume 13 - 2025 | https://doi.org/10.3389/feart.2025.1605436

Reverse time migration based on local Nyquist cross-correlation imaging condition with accurate first-arrival traveltimes correction

  • Sinopec Northwest China Petroleum Bureau, Urumqi, Xinjiang, China

Abstract

The selection of imaging conditions is one of the most critical factors determining the quality of reverse time migration (RTM) images. Among the widely used imaging conditions, the cross-correlation imaging condition (CCIC) consistently delivers high-resolution images. However, it is accompanied by substantial calculational costs and I/O tasks, particularly in 3D scenarios. In contrast, the excitation amplitude imaging condition (EAIC) offers advantages in computational efficiency, low storage requirements, and high precision. Nevertheless, it suffers from image distortion when dealing with multi-path propagation or strong reflection interfaces. The local Nyquist cross-correlation imaging condition (LNCIC) effectively combines the advantages of the two aforementioned imaging conditions. It uses the local wavefield near the time corresponding to the maximum amplitude at each grid point for imaging, and introduces the Nyquist sampling theorem to establish the search time step. This approach offers the benefit of high imaging quality while maintaining low storage cost. In this paper, we adopt an adaptive finite difference operator to solve the eikonal equation and calculate the accurate first-arrival traveltimes, thereby modify LNCIC and further enhancing the imaging accuracy. The effectiveness of the proposed method is demonstrated through numerical examples, including the Marmousi model, noise-resistance tests, and field data applications.

1 Introduction

As oil and gas exploration continues to advance, the focus of petroleum exploration has shifted towards complex geological formations and lithological exploration. As one of the most advanced techniques for using seismic data to describe subsurface structures and properties, reverse time migration (RTM) demands increasing imaging accuracy and computational efficiency in its algorithms. Since its inception in 1983 (; ; ), RTM has gained widespread acclaim among scholars due to its advantages, such as the no dip angle limitations, high imaging accuracy, and the ability to handle various events (; ; ).

The pre-stack RTM mainly involves three steps: first, the forward simulation to obtain the source wavefield; second, using the seismic records as boundary conditions to backpropagate and obtain the receiver wavefield; and finally, applying the imaging condition for imaging ().

Imaging conditions are a key step in influencing the quality of RTM images. The role of imaging condition is to select the optimal match between the wavefield predicted from the model space and the wavefield extrapolated from the data space to generate the image (). proposed that images are generated by multiplying two wavefields at each time step, known as the zero-lag cross-correlation imaging condition (ZLCCIC). However, this approach distorts the amplitude relationships of subsurface structures. To address this limitation, introduced the source-normalized cross-correlation imaging condition (SNCCIC), which normalizes the image using source illumination intensity. This method enhances imaging resolution while preserving reflectors, without introducing additional calculational costs. Although the SNCCIC can consistently deliver high-resolution images, its implementation requires storing the source wavefield at each time step to disk and subsequently reading it back during application. This approach imposes significant storage demands and severely impacts computational efficiency, particularly in 3D scenarios. summarized two strategies to reduce RTM storage requirements: the first one, the Nyquist method, which performs cross-correlation between source and receiver wavefields at Nyquist time intervals, and the second one, data compression techniques based on lossless compression algorithms. While these approaches can partially alleviate storage burdens, they still fall short of meeting practical application demands.

The wavefield reconstruction represents another class of approaches to address hardware storage limitations, primarily comprising two strategies. The first is the checkpoint method (), which selects some time steps as checkpoints. During wavefield propagation simulation, wavefield snapshots are stored in the computer’s memory or on disk at these checkpoints. These snapshots serve as initial conditions for reconstructing the wavefield through wave propagation between checkpoints. demonstrated that when the checkpoint intervals follow a binomial distribution, the computational efficiency is significantly improved. introduced a compute-for-storage strategy to optimize the checkpoint method, which requires storing only a small number of cached points to effectively reconstruct the wavefield. However, this approach necessitates repeated recursive wavefield reconstruction, with the number of recursions increasing exponentially as the number of time sampling points grows, resulting in a high recomputation rate. proposed a checkpoint technique for wavefield reconstruction based on interpolation principles. Under the condition of satisfying the Nyquist sampling theorem, they performed regular sampling of the wavefield between adjacent checkpoints, using the sampled wavefields as interpolation nodes. Polynomial interpolation algorithms were then applied to reconstruct the wavefield at any given time, optimizing the computational efficiency issues caused by the repeated iteration of the checkpoint technique. Another wavefield reconstruction approach is the boundary value method (; ), which stores boundary wavefields during forward modeling and uses these stored values to propagate the source wavefield backward from the final time step. While this method reduces storage requirements compared to full wavefield storage, its memory demand still increases with the order of the finite-difference operator. The introduction of random boundary conditions (; ; ) alleviates storage issues but introduces random noise, particularly in shallow regions near receivers. Furthermore, due to numerical instabilities associated with reverse-time reconstruction, this method cannot effectively handle attenuating media.

In practical applications, storing the wavefield on a hard drive or in memory is the simplest and most direct way to address the time consumption issue. Therefore, proposed the Excitation Amplitude Imaging Condition (EAIC). EAIC requires storing only the maximum amplitude of the source wavefield at each grid point as the excitation amplitude, along with the corresponding time as the excitation time. The imaging process involves dividing the receiver wavefield amplitude at the excitation time by the stored excitation amplitude. While EAIC offers inherent advantages in computational efficiency and wavefield storage requirements, it comes with certain trade-offs. Specifically, in scenarios involving multi-path wave propagation or complex media, the maximum amplitude may become distorted, leading to artifacts in the migrated images. By replacing excitation time with first-arrival traveltimes calculation in the improved EAIC, both imaging resolution enhancement and artifact suppression can be achieved (; ). The improved EAIC, while effective, remains limited by its utilization of only a subset of seismic data, resulting in inadequate suppression of background noise. To achieve both reduced storage requirements and high-quality imaging, proposed a local cross-correlation imaging condition (LCIC). This approach utilizes wavefields near the time of maximum amplitude at each grid point during forward modeling for imaging. Furthermore, they incorporated Nyquist sampling theory to develop a local Nyquist cross-correlation imaging condition (LNCIC), thereby further improving computational efficiency. However, the excitation time determination method in LNCIC shares the same limitation as EAIC - both rely on the time corresponding to the maximum amplitude of the source wavefield. This approach leads to highly disordered excitation time, which constitutes one of the primary reasons for EAIC’s poor imaging quality. To address this limitation, this study introduces a novel approach that employs an adaptive finite-difference method (AFDM) to accurately solve the eikonal equation for first-arrival traveltimes calculation (). Based on these traveltimes, we propose an enhanced LNCIC (ELNCIC) that utilizes the first-arrival traveltimes as the excitation time. This precision-based approach significantly improves the imaging accuracy of LNCIC.

The remainder of this paper is organized as follows: the section Ⅱ begins with a theoretical review of wavefield propagation in RTM and AFDM, followed by a detailed introduction to ELNCIC. In section Ⅲ, we validate the feasibility of ELNCIC using the Marmousi model, conduct noise-resistance tests, and demonstrate its practical applicability through field data examples. The section Ⅳ discusses the limitations of ELNCIC and an outlook on future developments, while section Ⅴ summarizes the paper.

2 Theory

2.1 The wavefield extrapolation in RTM

In this study, we only discuss acoustic wave propagation in isotropic media. The wave propagation is simulated using the second-order constant-density velocity-stress acoustic wave equation:where v is the velocity, x and z are the horizontal and vertical direction in 2D space, P (x, z) represents pressure, t is time. With the inclusion of the source term, the wavefield extrapolation in the forward direction can be expressed by Equation 2:where Ps(x, z) is source wavefield, represents Dirac function, xs and zs denote the source coordinates in the x- and z-directions, respectively, f(t) is source function. By incorporating seismic records into Equation 1, we obtain the expression for wavefield propagation in the reverse-time direction, as shown in Equation 3:where pᵣ represents the receiver wavefield, xᵣ and zᵣ denote the receiver coordinates in the x- and z-directions, respectively, and u indicates the seismic records.

We employ the finite-difference method to solve the wave equation, utilizing second-order temporal accuracy and sixth-order spatial accuracy, with convolutional perfectly matched layers (CPML) () as the absorbing boundary condition.

2.2 Review of adaptive finite-difference method

Under the high-frequency approximation, the seismic wave equation can be decomposed into an eikonal equation for traveltimes calculation and a transport equation for amplitude determination (; ). The traveltimes obtained by solving the eikonal equation plays a crucial role in seismic imaging (; ; ). The eikonal equation, a first-order nonlinear partial differential equation, can be derived from the Helmholtz equation in the high-frequency regime:where is the Hamilton operator, T(x) is traveltimes, S(x) is the slowness, Rn. For simplicity, we discretize the finite computational domain into grid points as illustrated, where the grid space in the x- and z-directions are denoted by dx and dz, respectively. To facilitate understanding, we have presented the distribution of traveltimes and slowness in the grid in Figure 1. Here, Tᵢ,ⱼ represents the numerical solution of first-arrival traveltimes at the grid vertices, while Sᵢ,ⱼ denotes the slowness at the grid cell centers. We adopt the AFDM to solve the eikonal equation. AFDM comprehensively accounts for all wave propagation modes by selecting the most accurate local operator—among plane-wave, spherical-wave, and refracted-wave operators—at each grid point to compute accurate first-arrival traveltimes.

FIGURE 1

2.2.1 Plane wave operator

In the finite-difference scheme, the partial derivatives of Equation 4 with respect to the x- and z-directions can be represented using the following central difference formulas:where Th and Tv are variables introduced to simplify subsequent complex expressions. By incorporating Equations 5, 6 into Equation 4, followed by straightforward mathematical derivation, the first-arrival traveltimes for the plane wave can be determined as follows:

2.2.2 Spherical wave operator

To overcome the singularity in simulating first-arrival traveltimes near the source, the solution for T(x) can be expressed as the product of a known solution and a decomposition factor. The traveltimes and slowness model are decomposed according to Equation 8:

If both T0(x) and S0(x) are known solutions, Equation 4 can be transformed to numerical evaluation of the factor , as given by Equation 9:Equation 10 defines the boundary conditions:we can calculate a numerical solution for the factor τ(x), τh and τv can be obtained by replacing T with τ in Equation 6:

2.2.3 Refracted wave operator

The refracted wave operator has been developed to accurately simulate the propagation of refracted or direct waves:

The AFDM calculates first-arrival traveltimes at each grid point using the plane wave operator from Equation 7, the spherical wave operator from Equation 11, and the refracted wave operator from Equation 12. It selects the minimum value among these three as the candidate traveltimes Tcan. If Tcan is smaller than the previous one, the method updates the traveltimes accordingly. This process continues until the traveltimes converge to stable values.

2.2.4 Accuracy analysis of adaptive finite-difference method

To validate the accuracy of first-arrival traveltimes computed by the AFDM, we compare AFDM with Godunov upwind difference method () and the hybrid method (). Since the analytical solution for traveltimes from a source in a homogeneous model can be easily derived, we construct a 60 × 60 grid homogeneous slowness model with a slowness of 0.0004 s/m. The relative errors between the numerical solutions and the analytical solution for the three methods are calculated using the Equation 13:where Re is relative error, Snum and Sana represent the numerical solution and the analytical solution, respectively.

Figure 2 presents the relative errors between the numerical solutions obtained by the three methods and the analytical solution. Due to source singularity, Figure 2A shows noticeable relative errors near the source for Godunov upwind difference method. In Figure 2B, the hybrid method exhibits smaller relative errors near the source owing to its windowing scheme; however, it fails to adequately handle wave propagation outside the window, resulting in significantly increased errors in those regions. In contrast, Figure 2C demonstrates that AFDM achieves substantially smaller maximum relative errors compared to the other two methods, confirming its accuracy in computing first-arrival traveltimes. Table 1 presents a comparative analysis of computational time and maximum relative errors for three different eikonal equation solvers. While AFDM requires slightly longer computation time than the other two methods - owing to its utilization of three different local operators for eikonal equation solutions - it achieves significantly smaller maximum relative errors. We contend that this marginal increase in computational overhead is well justified by the substantial improvement in first-arrival traveltimes accuracy.

FIGURE 2

TABLE 1

MethodsGodunov upwind difference methodHybrid methodAFDM
Elapsed time (s)6.1e-35.6e-38.09e-3
Maximum relative error (%)0.214.4e-31.5e-4

Comparison of elapsed time and maximum relative error for three different eikonal wave equation solvers.

2.2.5 The ELNCIC based on adaptive finite-difference method

Although the first-arrival traveltimes computed using the AFDM are sufficiently accurate, the eikonal equation calculation inherently ignores the wavelet duration. Additionally, considering various sources of error in practical applications, this study adopts the following reference first-arrival traveltimes Tre is determined by Equation 14:where TAFDM is first-arrival traveltimes calculated by AFDM, fm represents dominant frequency of source wavelet, a is an error factor introduced to account for practical considerations ().

In the LNCIC framework (), each grid point maintains a time window of length 2L+1 to store local wavefield. During forward and backward propagation, we preserve the source and receiver wavefield within their respective time windows, which are subsequently applied to imaging using Equation 15:where ILNCIC represents the migration image of the LNCIC, sj denotes the jth shot, Tmax refers to the excitation time, and ti is the ith time step, ts is the search step length proposed based on the Nyquist sampling theorem in Equation 16:

where fmax and fmin represent the maximum and minimum frequency of the seismic records, respectively. To demonstrate that the excitation time used in LNCIC affects the quality of the migration image, we use the Marmousi velocity model (Figure 3) and set the source position at 2062.5 m to obtain the excitation time and single-shot migration image for analysis. Figure 4A shows the excitation time used in LNCIC, and Figure 4B presents the single-shot migration image obtained by cross-correlating the wavefield corresponding to this excitation time. It is evident in the figure that there is noticeable distortion in the migration image. We have also used arrows to highlight the locations of distortion in the migration image, which directly correspond to the disturbed positions of the excitation time. Therefore, using a more continuous and accurate first-arrival traveltimes as the excitation time will effectively improve the accuracy of the migration image. At the same time, to demonstrate the feasibility of replacing the excitation time with the first-arrival traveltimes calculated by AFDM, we set the source at the same location to obtain the first-arrival traveltimes calculated by AFDM (Figure 5A) and the single-shot migration image calculated using this time (Figure 5B). We can observe that, due to the stability and continuity of the calculated first-arrival traveltimes, there is no distortion in the migration image in Figure 5B.

FIGURE 3

FIGURE 4

FIGURE 5

By incorporating Tre into LNCIC, we can obtain the following expression for the ELNCIC from Equation 17:where IELNCIC represents migration image of the ELNCIC. The stacked migration image is obtained by stacking the migration images from each shot.

3 Numerical examples

In the numerical examples, we compare the performance of ELNCIC with three methods: the traditional EAIC, the stable EAIC (SEAIC) () and SNCCIC. These comparisons are conducted using the Marmousi model, noise-contaminated data, and field data to evaluate imaging quality and robustness.

3.1 Marmousi model example

The Marmousi model has a grid size of 663 × 234 with a spatial sampling interval of 6.25 m (Figure 3). A total of 100 shots are evenly distributed along the surface, with a shot spacing of 37.5 m. A Ricker wavelet with a dominant frequency of 20 Hz is used as the source, and the recording duration is 3 s with a temporal sampling interval of 0.5 m. For the selection of the time window size, we use the period of the source wavelet, tsource, as a measure. Through our tests, we found that choosing a time window size of 3.0tsource provides high-quality imaging while maintaining relatively small computational time and disk storage requirements. The numerical examples in this paper all use this window size.

In Figure 6, we present the migration results obtained using four different imaging conditions. The images generated by EAIC (Figure 6a) exhibit significant distortion in the red rectangular region. This artifact arises from the inherent limitations of EAIC, where the maximum amplitude values in these areas may become distorted because the stored excitation amplitudes do not correspond to first-arrival waves but rather to other wave types, such as reflected or refracted waves. The resulting amplitude inconsistencies introduce severe interference during multi-shot stacking. Figure 6b demonstrates that SEAIC, which employs precise first-arrival traveltimes as excitation times, significantly reduces the distortion within the red rectangle. The migration image generated by ELNCIC (Figure 6c) not only completely resolves distortion in complex structural areas but also effectively suppresses background noise, demonstrating comparable quality to the SNCCIC result (Figure 6d). However, ELNCIC achieves this while significantly reducing storage requirements and improving computational efficiency, as it only utilizes the local source wavefield and receiver wavefield, unlike SNCCIC which requires full wavefield storage.

FIGURE 6

To facilitate clearer comparison and observation, we provide an enlarged view of the red rectangular region in Figure 7. Additionally, we extract two waveform traces from the migration images at x = 312.5 m and x = 1562.5 m (Figure 8). These traces provide a more intuitive comparison of the four imaging conditions. The ELNCIC results show excellent agreement with SNCCIC, while black arrows highlight the distorted regions in the EAIC and SEAIC results.

FIGURE 7

FIGURE 8

Table 2 shows the disk storage space and computation time required for the four imaging conditions. From the data in the table, it is evident that ELNCIC requires only 2.23% of the storage and 21.99% of the computation time compared to SNCCIC. Although the storage requirements and computation time for ELNCIC are slightly higher than those for EAIC and SEAIC, the resulting image quality more than justifies the increase.

TABLE 2

Imaging conditionsStorage costs (Mb)Elapsed time (s)
EAIC1.11101.55
SEAIC59.90115.12
ELNCIC148.56135.29
SNCCIC6645.76615.10

Comparison of storage cost and elapsed time for the four imaging conditions.

3.2 Noise resistance analysis

In practice, seismic data are often severely contaminated by noise, and denoising during preprocessing may result in the loss of valuable signals. To preserve the original seismic data as much as possible, we add random Gaussian noise to the seismic records to reduce the signal-to-noise ratio (SNR). The Equation 18 gives the expression for SNR:where, Psignal and Pnoise are represent the signal power and noise power, respectively. We set the SNR of the seismic records to 5 dB. Figure 9 compares the original seismic records and seismic records with SNR = 5 dB, revealing that some valuable signals in the seismic records are severely contaminated.

FIGURE 9

Figure 10 demonstrates the imaging performance of the four imaging conditions on noise-contaminated seismic data. EAIC exhibits poor noise resistance, causing severe contamination in the migration images, particularly in deeper layers. Although SEAIC has been improved by utilizing the highest SNR portions of the seismic data for imaging, it still fails to adequately suppress background noise. In contrast, both ELNCIC and SNCCIC significantly outperform the former two methods. ELNCIC, which not only effectively suppresses noise and enhances image quality but also substantially reduces storage requirements and improves computational efficiency compared to SNCCIC, which uses wavefields from all time steps.

FIGURE 10

Similarly, we extracted waveform traces at x = 312.5 m and x = 1562.5 m (Figures 11a,b), using the black curve representing SNCCIC as the reference. Under noise contamination, both EAIC (green curve) and SEAIC (blue curve) exhibit significantly distorted waveforms that fail to match the reference, whereas ELNCCIC (red curve) maintains good agreement with SNCCIC, demonstrating its superior noise resistance capability.

FIGURE 11

3.3 Field data example

Next, we evaluate the performance of ELNCIC on field data. The velocity profile, shown in Figure 12, spans a width of 15.61 km and a depth of 5.5 km. The computational domain is discretized into a 1561 × 550 grid with a spatial sampling interval of 10 m. A total of 88 shots are unevenly distributed along the surface, with each shot recorded for 6 s at a temporal sampling interval of 0.5 m.

FIGURE 12

Figure 13 demonstrates the application of AFDM in field data. The calculated first-arrival traveltimes do not exhibit any disturbances, indicating that this method is also effective in field data. Figure 14 demonstrates the performance of the four imaging conditions on field data. The regions within the red and blue rectangles provide a clearer comparison of the differences among the four methods. Due to their reliance on wavefields only at excitation times, EAIC and SEAIC produce migration images with weaker event amplitudes, particularly below 4 km depth, where their performance is unsatisfactory. In contrast, ELNCIC, as theoretically predicted, generates images with continuous and clear stratigraphic features, achieving imaging quality comparable to SNCCIC. For a more detailed observation, enlarged views of the red rectangular regions in Figure 14 are presented in Figure 15.

FIGURE 13

FIGURE 14

FIGURE 15

4 Discussion

The ELNCIC proposed in this study, similar to EAIC and SEAIC, belongs to a specific category of cross-correlation-based imaging conditions. While EAIC and SEAIC utilize wavefield information only at excitation times for cross-correlation imaging, ELNCIC employs locally effective wavefield information. Theoretically, richer wavefield information used for cross-correlation leads to higher-quality migration images, but this comes at the cost of substantial storage requirements and high computational expenses. EAIC and SEAIC, which utilize limited wavefield information, achieve computational efficiency but offer virtually no noise resistance. Therefore, ELNCIC represents an effective compromise, delivering high-quality imaging and strong noise resistance while maintaining reasonable computational efficiency.

When addressing multi-path issues, ELNCIC can adapt by appropriately increasing the time window length, trading some computational efficiency for improved performance in complex wave propagation scenarios. While this study focuses on isotropic media, future work will extend the application to more complex media, including attenuative and anisotropic media.

5 Conclusion

The ELNCIC uses the first-arrival traveltimes, calculated by AFDM, as the excitation time. A time window is set around the excitation time, and the Nyquist sampling theorem is applied to establish the search time step. By storing the local source and receiver wavefields for imaging, it achieves low storage requirements while maintaining high-quality imaging. Through numerical experiments on the Marmousi model, noise-contaminated data, and field data, we have demonstrated ELNCIC’s superior imaging capabilities in complex geological settings and its robustness against noise.

Statements

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

GL: Writing – original draft, Data curation, Methodology, Writing – review and editing. LZ: Funding acquisition, Writing – review and editing. LH: Writing – review and editing. YW: Writing – review and editing. ZQ: Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work is supported by the Study on identification, description and evaluation technology of ultra-deep and super-deep carbonate rock traps (P24009) and Integrated research and application of ultra-deep and super-deep drilling geological engineering (P24136).

Acknowledgments

We would like to thank the reviewers and editors for their valuable suggestions for the paper.

Conflict of interest

Authors GL, LZ, LH, YW, and ZQ were employed by Sinopec Northwest China Petroleum Bureau.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Summary

Keywords

reverse time migration, adaptive finite- different method, eikonal equation, first-arrival traveltimes, imaging conditions

Citation

LiJun G, ZongJie L, HaiYing L, Wei Y and Qing Z (2025) Reverse time migration based on local Nyquist cross-correlation imaging condition with accurate first-arrival traveltimes correction. Front. Earth Sci. 13:1605436. doi: 10.3389/feart.2025.1605436

Received

03 April 2025

Accepted

23 April 2025

Published

12 May 2025

Volume

13 - 2025

Edited by

Jidong Yang, China University of Petroleum (East China), China

Reviewed by

Hanming Chen, China University of Petroleum, China

Yong Hu, China University of Mining and Technology, China

Yufeng Wang, China University of Geosciences Wuhan, China

Updates

Copyright

*Correspondence: Gao LiJun,

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics