The advancement of the petroleum industry has heightened the demand for high-precision imaging techniques in seismic exploration. Seismic migration, like taking a photograph of underground structures, resolves geometric distortions and noise interference in raw seismic data. As an indispensable technique in seismic exploration, its essence is to image through the convergence points of the forward wavefield and the reverse wavefield, providing reliable foundations for geological interpretation. It primarily includes Kirchhoff Migration (KM), Reverse Time Migration (RTM) and Gaussian Beam Migration (GBM). KM, characterized by its computational efficiency and operational flexibility, serves as a widely adopted tool in seismic imaging. Particularly effective for rapidly delineating medium-complexity geological structures, this method fundamentally employs ray-theoretical approximations and integral stacking methods, constituting a speed-for-accuracy compromise strategy. In contrast, RTM utilizes two-way wave equation solutions to achieve high-fidelity subsurface imaging. By synthesizing forward source wavefield with back receiver wavefield through imaging conditions, it demonstrates superior performance in resolving intricate geological features. As a high-frequency asymptotic solution to the wave equation, GBM satisfies both regularity and concentration conditions, achieving imaging accuracy comparable to RTM. Developed from ray theory, GBM combines high flexibility with computational efficiency. Additionally, it adheres to the regularity condition, avoiding amplitude singularities and thus eliminating imaging issues associated with caustics and shadow regions. Due to these unique advantages, GBM has been the subject of extensive research and is increasingly being implemented in practical applications.

The Gaussian beam method has undergone systematic theoretical evolution in seismic wavefield simulation and imaging applications since its inception. The foundational work was established by Červený et al. (1982), who pioneered the application of Gaussian beam superposition theory to seismic wavefield modeling, laying the theoretical groundwork for subsequent developments in Gaussian beam forward modeling and migration. Building upon this framework, Hill (1990) and Hill (2001) made significant advancements by developing the first frequency-domain GBM (FGBM) imaging system and introducing the concept of pre-stack depth migration. His work demonstrated that it maintains the computational advantages of ray-based migration while exhibiting superior adaptability to complex wave phenomena, including wavefront dispersion and shadow zone effects. To address inherent limitations in pre-stack FGBM implementations, Gray (2005) proposed an optimized common shot domain FGBM approach. This innovation not only enhanced azimuthal data processing capabilities but also provided more effective solutions for near-surface seismic challenges. Subsequent technical refinements emerged through the work of Hu and Stoffa (2009), who derived a slowness-driven FGBM formulation incorporating horizontal surface slowness parameters. Their methodology achieved natural integration of Fresnel weighting with beam superposition, effectively suppressing artifacts caused by incomplete superposition and thereby improving migration fidelity. Yang and Zhu (2018a) developed a data-driven optimization strategy to mitigate offset artifacts prevalent in low signal-to-noise ratio datasets, significantly enhancing interpretability. Through systematic analysis of crust-mantle structures, Han et al. (2022) quantitatively validated FGBM as a high-precision pre-stack depth migration method for deep reflection seismic data, establishing critical theoretical benchmarks for practical implementations. Recently, Yu et al. (2022) extended the method’s applicability through vertical seismic profile implementations, introducing a free-surface multiple FGBM technique that effectively compensates for shallow structural imaging energy loss in far-offset acquisitions.

The application of GBM in complex topographic environments has undergone significant methodological advancements to address surface-related imaging challenges. Initial breakthroughs emerged from the work of Yue et al. (2010) and Yue et al. (2012), who developed the amplitude-preserved FGBM formulation for rugged terrains through fundamental modifications to the tilt-stacking operator. By explicitly incorporating surface elevation parameters into the plane-wave synthesis algorithm, their method achieved enhanced accuracy in local wavefield reconstruction under complex surface geometries. Subsequent field applications demonstrated the method’s practical value in mountainous regions. Cao et al. (2013) successfully implemented this amplitude-preserved FGBM algorithm in piedmont zone environments, providing the systematic validation of its adaptability to fields with complex topographic structures. Building upon these topographic adaptations, Huang et al. (2015a) introduced a FGBM framework specifically designed for dual-complexity scenarios involving both irregular surface geometries and heterogeneous subsurface media. They maintained computational efficiency while handling multi-scale wave propagation effects. Recently, Han et al. (2020) extended the FGBM framework of Huang et al. (2015a) to anisotropic elastic media through rigorous incorporation of Thomsen anisotropy parameters. Their anisotropic elastic FGBM formulation not only preserved the method’s inherent advantages in handling surface complexity but also enabled more physically accurate imaging in stratified media with directional velocity variations.

The extension of GBM to complex media has driven critical methodological innovations to address wave propagation complexity in realistic subsurface environments. Initial theoretical breakthroughs emerged from Alkhalifah’s (1995) FGBM formulation for transversely isotropic (TI) media, which introduced anisotropic ray-tracing modifications but required precise estimation of Thomsen parameters. To overcome this limitation, Zhu et al. (2007) rederived the ray-tracing system using phase velocity parametrization, establishing a generalized TI media framework with reduced parameter sensitivity. Building upon these foundations, Han et al. (2014) systematically quantified P/SV-wave propagation in 2D TI media through kinematic/dynamic ray tracing, coupled with PS-wave common-shot imaging conditions, thereby enabling converted-wave pre-stack depth migration in anisotropic media. Bai et al. (2016) derived viscoacoustic Green’s functions via Gaussian beam superposition, incorporating Q-compensation mechanisms to address the attenuation effects. Concurrently, Li and Mao (2016) resolved 3D polarity reversals in elastic wavefield by developing cross-correlation imaging conditions integrated with beam-based wavefield extrapolation. Subsequent refinements by Li et al. (2018) introduced converted-wave ray tracing in vertical transversely isotropic (VTI) media, while Xu and Mao (2018) enhanced computational efficiency through multilayer anisotropic ray tracing with complex initial conditions, explicitly accounting for multipath effects through rigorous model validation. Some other researches have focused on multi-physics integration and practical optimization. Yang et al. (2018b) achieved simultaneous P/S-wavefield extrapolation via Kirchhoff-Helmholtz decomposition, paired with an enhanced dot-product imaging condition to improve PS-wave resolution. Qin (2020) further streamlined anisotropic implementations through GBM formulations with optimized ray-tracing coefficients. Recently, Shi et al. (2023) addressed sparse ocean-bottom node imaging challenges by developing a 4-component elastic FGBM framework based on elastic reciprocity theory, demonstrating robust performance in low-coverage acquisition scenarios.

The continuous refinement of GBM has driven systematic improvements in imaging precision through beam optimization strategies. Initial progress was marked by Nowack’s (2008) and Nowack’s (2011) dynamic focused beam migration framework, which introduced spatial control of beam geometry by concentrating the narrowest beam segment within predefined target zones, thereby establishing the foundational principle of beam focusing for resolution enhancement. Building upon this foundation, Yang et al. (2015a) developed an amplitude-preserved adaptive focusing beam method that synergistically integrated multiple focusing operators with velocity-dependent beam waist adjustments, which is like a camera focusing, thus improving the imaging accuracy of the target area. Concurrently, Huang et al. (2015b) proposed a paradigm shift through Fresnel beam migration, constraining beam geometry via Fresnel zone width limitations rather than conventional Gaussian beam. Their methodology achieved superior wavefield reconstruction accuracy, particularly in complex near-surface environments, while expanding applicability to cross-well seismic configurations through optimized beam superposition criteria (Yang et al., 2022). Recent advancements have extended these beam optimizations to complex media. Liu et al. (2022a) extended the focused beam migration algorithm to anisotropic media, and used anisotropic ray tracing equations with optimized coefficients to obtain the beam information, and proposed an anisotropic dynamic focused beam migration method. Similarly, focused beam migration is also applied to viscoacoustic media (Zhang et al., 2024).

In order to further optimize the results of migration imaging, many scholars have carried out the Least Square GBM (LSGBM) method. Hu et al. (2016) first proposed the concept, compared with the traditional GBM method, the resolution and amplitude fidelity of the migration results can be further improved through continuous iterative optimization. Based on the Born approximation theory of isotropic acoustic wave equation, Yang et al. (2018c) derived a time-domain linear Gaussian beam modeling operator, and proposed a time-domain LSGBM method, which greatly reduced the computational cost. Since then, the LSGBM method has gradually considered the influence of elasticity and viscosity of the medium (Yue et al., 2019; Yue et al., 2021; Mao et al., 2022).

Space-time-domain GBM (STGBM) is a higher precision seismic imaging method well developed in recent years compared to FGBM. STGBM utilizes Gaussian beams to model the propagation of seismic wave in the space-time domain, enabling the construction of both forward and reverse seismic wavefield for migration imaging. This method primarily encompasses two techniques: Gaussian Envelope Migration and Gaussian Beam Reverse Time Migration (GBRTM). The foundational work in this field was conducted by Žáček (2006), who first realized Gaussian Envelope Migration by deriving a series of space-time-domain wavefield formulas based on Gaussian envelope superposition. Building on this, Popov et al. (2010) developed a GBRTM method using Kirchhoff integral-based reverse time extrapolation, laying the theoretical groundwork for STGBM. A significant advancement was made by Yang et al. (2015b), who implemented wavefield back-propagation using Gaussian beams and employed Gabor transform for sparse decomposition of seismic wavefield. This STGBM approach demonstrated superior performance with reduced migration artifacts compared to FGBM. Further developments in STGBM have addressed various complexities in seismic imaging. Huang et al. (2017) extended the method to account for elastic media characteristics by applying Gaussian beam superposition to RTM imaging. Lv et al. (2019) enhanced the method’s performance in structures with sharp lateral velocity variations through the introduction of dynamic focusing theory. The method’s applicability to elastic media was further improved by Hu et al. (2020), who incorporated wavefield separation of pure P and S waves in elastic STGBM. Recent advancements have focused on computational efficiency and practical applications. Zhang et al. (2022a) and Zhang et al. (2022b) developed a fast STGBM method using dominant frequency approximation and ray optimization strategies, significantly improving computational efficiency while maintaining accuracy in VTI media. Recently, Liu et al. (2022b) proposed an optimized STGBM method for TI media through enhanced ray tracing techniques.

When the field is dealing with the deep exploration data in the large work area, especially the 3D massive data, it is optimal to use GBM. Therefore, it is very important for engineers and researchers to study GBM methods and applications. In this paper, we will first introduce in detail the FGBM and STGBM methods in acoustic media, so as to enable readers to have a good understanding of the principles and implementation methods of the two GBM methods, and then demonstrate the imaging differences between the two types of methods using simulated data. Furthermore, we will use two STGBM methods respectively in acoustic and VTI media to image the actual oil field data, aiming to help engineers and researchers better understand the application of GBM and the necessity of taking into account the real mechanical properties of the media. Finally, we will present the calculation times of different methods under three classic scenario conditions, and based on these results, we will provide specific suggestions for selecting different imaging methods.

In this chapter, we mainly present the implementation of the FGBM and some important formulas. The detailed derivations are from the paper (Červený et al., 1982). In 2D ray-centered coordinate system (Figure 1), the travel time along a ray path can be expressed as: τ s = ∫ 0 s 1 v s d s (1) where τ s is the travel time along the central ray, s is the arc length along the central ray, and v s is the velocity of the seismic wave.

In acoustic media, consider the following scalar wave equation in 2D Cartesian coordinate system: Δ U x , z , ω + ω 2 V x , z 2 U x , z , ω = 0 (2) where x , z represents the coordinates of a 2D space point, ω is the angular frequency of the seismic wave, U x , z , ω represents the displacement of the seismic wavefield related to space and frequency, V x , z is the velocity at x , z , and Δ is Laplacian operator, which can be written in the following form within 2D Cartesian coordinate system: Δ = ∇ 2 = ∂ ∂ x 2 + ∂ ∂ z 2 (3)

Numerous methods exist for solving wave equations, including ray tracing, finite difference method, finite element method, and spectral element method. Among these, ray tracing offers an efficient kinematic and dynamic description of the wavefield, circumventing the computational expense of full wave equation solvers while effectively capturing wave propagation in complex media. Building upon the ray-centered coordinate system, Červený et al. (1982) formulated Equation 2 and derived the ray tracing system of equations and the high-frequency asymptotic solution to the wavefield through high-frequency approximations and variable substitutions, establishing the theoretical framework for Gaussian beam forward modeling. Extending this work, Hill (1990) developed a migration theory based on Gaussian beam superposition, employing different notation to reexpress the same forward modeling equations as those presented by Červený et al. (1982).

In 2D problems, the wave equation for ray tracing is given by: d Q s d s = v s P s d P s d s = − 1 v 2 s ∂ V s , n ∂ n 2 Q s (4) where v s = V x s , z s , P s and Q s are complex solutions of the ray tracing equations. Specifically, P s characterizes the dynamic properties of the beam, including amplitude attenuation and phase variation, while Q s describes the geometric attributes such as beam width and curvature.

The solution to Equation 4, corresponding to the Gaussian beam expression, is given by: u GB s , n , ω = v s Q s exp i ω τ s + 1 2 i ω P s Q s n 2 (5)

Figure 2A shows the ray field in a simple two-layer model. The ray field is calculated using Equation 4, and the two-layer model is discretized into a 201 × 201 grid, with a sampling interval of 10 m in both horizontal and vertical directions. It is easy to see that the directions of the rays change when they hit the interface. Figure 2B is the Gaussian beam of the two-layer model at the reference frequency of 20 Hz.

Meanwhile, the seismic receiver wavefield U r , ω can be expressed via the integral relationship between the seismic source wavefield U r ′ , ω and the complex conjugate of the Green’s function G * r , r ′ , ω : U r , ω = − 1 2 π ∬ d x ′ ∂ G * r , r ′ , ω ∂ z ′ U r ′ , ω (6) where r = x , z and r ′ = x ′ , z ′ represent the coordinates of the receiver and shot points, respectively. G is the Green’s function and G * represents its complex conjugate. The Green’s function can be expressed as the superposition of a series of Gaussian beams (Hill, 2001): G r , r ′ ; ω ≈ i 2 π ω ∬ d p x ′ p z ′ u GB r ; r ′ , p ′ ; ω (7) where p ′ = p x ′ , p z ′ is the slowness of the shot point.

According to the imaging principle, the formula of the FGBM can be derived from the cross-correlation between the seismic receiver wavefield U r , ω and the seismic source wavefield U r ′ , ω : I FGBM = ∫ d ω U r , ω U r ′ , ω (8)

STGBM employs a hybrid strategy integrating wave equation approximations with ray-tracing techniques to enhance imaging accuracy at the expense of increased computational overhead. This methodology establishes an accuracy-complexity equilibrium, positioning itself as a balance between FGBM and RTM. In the space-time domain, the expression of the Gaussian beam can be written in the following form (Popov et al., 2010): U x 0 , t ; x s = − i 4 π ∫ d ω ε G B s 0 v s v s 0 Q s exp i ω − t − τ + 1 2 P s Q s n 2 (9) where x 0 = x 0 , z 0 and x s = x s , 0 represent the locations of the imaging and shot points, respectively, and ε G B is the initial Gaussian beam parameter.

According to the Gaussian beam forward modeling theory (Červený et al., 1982), the forward wavefield can be represented by a series of Gaussian beams emitted in different angles: W 1 = − i 4 π ∫ 0 2 π d φ ∫ d ω ε G B s 0 v s v s 0 Q s exp i ω − t − τ + 1 2 P s Q s n 2 (10) where φ is the angle at which the Gaussian beam emerges from the surface.

The back propagation process of the recorded wavefield from the receiver point to the imaging point can be realized by the Kirchhoff integral (Popov et al., 2010). In other words, the up-going ray tracing strategy is used to construct the reverse wavefield, as shown in Figure 3: W 2 = − 2 ∫ 0 T d t ∫ d x r P U x r , t ∂ ∂ z G x r , t − t 0 ; x 0 (11) where P U x r , t is the recorded wavefield, which is used as the input in the migration algorithm; G x r , t ; x 0 , t 0 is Green’s function, which can be expressed in the following form: G ⁡ x r , t − t 0 ; x 0 = − i 4 π ∫ 0 2 π d φ ∫ d ω ε G B s 0 v s v s 0 Q s ⁢ exp ⁡ i ω ⁡ − t − t 0 + 1 2 ⁢ P s Q s ⁢ n 2 (12)

In the migration algorithm, the imaging condition is very important, which can directly affect the final image quality. In the time window t ∈ t 1 , t 2 , when the wavelet phases of the forward wavefield W 1 and the reverse wavefield W 2 are the same, the cross-correlation of the two reaches the maximum value, the formula of the STGBM method in acoustic media can be obtained: I STGBM = ∫ t 1 t 2 d t W 1 W 2 (13)

We construct a multi-layer model to verify the above two GBM methods, as shown in Figure 4. Figure 4A is the realistic velocity model, and Figure 4B is its smoothed version. In the forward modeling process, an observation system with middle excitation and two sides receiving is adopted. A total of 201 shots are set up with a spacing of 30 m, the position coordinate of the first shot is (2 km, 0), and the position coordinate of the last shot is (8 km, 0). The dominant frequency of the wavelets is 20 Hz. Figure 5 shows the results of the first 1,200 traces of the shot records.

Figure 6 shows the imaging results using the above two GBM methods. Both methods can produce relatively clear migration results and image the karst cave structures in the middle layer, which verifies the effectiveness of the above two GBM methods. As shown by the red arrows in Figure 6, the FGBM method (Hill, 1990; Hill, 2001; Gray, 2005) produces some migration artifacts near the shallow layer, while the STGBM method (Yang et al., 2015b; Lv et al., 2019; Zhang et al., 2022a) is clearer and sharper. We extract the longitudinal single-trace results of the two in the shallow layer when the transverse distance is 4 km, it can be obviously found that the STGBM results have higher resolution in the shallow part (Figure 6C). This is because the STGBM computes more beam information at the same time, which improves the imaging accuracy of the migration results at a higher computational cost.

In seismic exploration, subsurface media are often complex. The VTI model is representative for most sedimentary basins with relatively low parameter complexity. In terms of elastic parameters, there are three parameters for acoustic VTI media: Vp (vertical P-wave velocity), ε (the difference between horizontal and vertical P-wave velocities), and δ (the curvature of the P-wave velocity variation with the angle). For elastic VTI media, there are five parameters in total, with the addition of Vs (vertical S-wave velocity) and γ (the difference between horizontal and vertical S-wave velocities). Due to these characteristics of the VTI model, we will mainly introduce the principle and implementation of the GBM in VTI media, and compare the imaging differences between the VTI STGBM and acoustic STGBM.

In VTI media, the implementation of GBM mainly modifies the ray tracing parameters. Hanyga (1986) presented the dynamic ray-tracing equations in VTI media: d Q s d τ = W Q s + V P s d P s d τ = − H Q s − V P s (14) where W, V and H are coefficient matrices related to the elastic properties.

The Gaussian beam propagation equation, reverse wavefield construction and migration are similar to the GBM methods in acoustic media.

Next, we present results using the actual data from an oil field to verify the STGBM method in VTI media. Figure 7A shows the velocity field, and the anisotropic parameters ε and δ are shown in Figures 7B, C, respectively. In the process of seismic data acquisition, the observation mode of unilateral reception is adopted. A total of 204 shots are set up, each shot has 60 traces, the position coordinate of the first shot is (3 km, 0), and the position coordinate of the last shot is (13.2 km, 0); both shot spacing and trace spacing are 50 m. The first 300 traces of the shot records for the actual data are shown in Figure 8.

Figure 9 shows two imaging results obtained by the acoustic STGBM (Yang et al., 2015b; Lv et al., 2019; Zhang et al., 2022a) and the VTI STGBM (Zhang et al., 2022b; Liu et al., 2022b), and Figures 10, 11 are the magnified views of the red and blue boxes in Figure 9, respectively. Compared with the blue ellipse, it can be seen that the migration results generated by the GBM method in VTI media exhibit relatively clear fault characteristics, which provides a more favourable condition for the subsequent oil reservoir interpretation. This is because the VTI GBM method considers the influence of the anisotropy, and more accurately calculates the travel time and amplitude of the seismic beam, so that the reflected wave of the offset result is correctly located.

We systematically introduce GBM methods from the perspectives of media properties and computational domains, supplemented by examples utilized three GBM methods. A multi-layer model is employed to show the effectiveness of the GBM, demonstrating that both FGBM and STGBM can effectively image complex structures; however, the STGBM exhibits reduced noise in shallow regions compared to the FGBM. Then, we separately apply acoustic and VTI STGBM to process a field dataset; the results indicate that the VTI GBM method provides superior fault image clarity compared to the acoustic GBM method.

In this study, we demonstrate the time consumption of FGBM and STGBM through the calculation of a multi-layer model. The FGBM demonstrates computational efficiency, with processing times per shot on a single-core CPU on the order of 10^–2 ∼ 10^–1 h for an array with dimensions of 401 (traces) × 1,001 (time sampling points); in contrast, the STGBM requires longer computational durations, typically on the order of 10^–1 ∼ 10⁰ h. The oilfield sector currently encounters significant challenges in processing ultra-large-scale work area datasets. As shown in Table 1, we present the order of magnitude of computation time per shot for the FGBM, STGBM, and RTM methods under different classic scenarios on a single-core CPU. In these scenarios, RTM incurs prohibitively high computational costs, where GBM presents itself as the optimal solution, striking an effective balance between computational efficiency and imaging accuracy. During the process of selecting the GBM, it is advisable to prioritize the application of the FGBM for its computational efficiency. However, in scenarios where substantial computational resources are accessible, including the utilization of supercomputing facilities, the adoption of the STGBM is recommended to achieve enhanced accuracy. When addressing diverse and complex work area data, it is essential to approach the problem by considering medium properties, such as anisotropy, elasticity, and viscosity, to progressively enhance the imaging quality of migration results (Alkhalifah, 1995; Bai et al., 2016; Shi et al., 2023). Furthermore, in cases where sufficient computational resources are available, the LSGBM (Hu et al., 2016; Yang et al., 2018c; Yue et al., 2019, 2021; Mao et al., 2022) can be employed to achieve higher resolution and improved fidelity in the migration results.

The GBM represents a versatile imaging tool that not only preserves the advantages of traditional ray-based migration but also achieves imaging accuracy comparable to wave-equation migration. As such, it serves as a highly effective approach for seismic data processing. We propose that its future development will primarily focus on the following directions:

1) Research on GBM methods in more complex media for deep structure areas;