Kirchhoff Migration Research Articles

Research on Airborne Ground-Penetrating Radar Imaging Technology in Complex Terrain

The integration of ground-penetrating radar (GPR) with unmanned aerial vehicles (UAVs) enables efficient non-contact detection, performing exceptionally well in complex terrains and extreme environments. However, challenges in data processing and interpretation remain significant obstacles to fully utilizing this technology. To mitigate the effects of numerical dispersion, this paper develops a high-order finite-difference time-domain (FDTD(2,4)) three-dimensional code suitable for airborne GPR numerical simulations. The simulation results are compared with traditional FDTD methods, validating the accuracy of the proposed approach. Additionally, a Kirchhoff migration algorithm that considers the influence of the air layer is developed for airborne GPR. Different processing strategies are applied to flat and undulating terrain models, significantly improving the identification of shallowly buried targets. Particularly under undulating terrain conditions, the energy ratio method is introduced, effectively suppressing the interference of surface reflections caused by terrain variations. This innovative approach offers a new technical pathway for efficient GPR data processing in complex terrains. The study provides new insights and methods for the practical application of airborne GPR.

Application of Kirchhoff Migration from Two-Dimensional Fresnel Dataset by Converting Unavailable Data into a Constant

Application of Kirchhoff Migration from Two-Dimensional Fresnel Dataset by Converting Unavailable Data into a Constant

Efficient pure qP-wave modeling and reverse time migration in tilted transversely isotropic media calculated by a finite-difference approach

Subsurface anisotropy is commonly induced by shale layers, aligned cracks, and fine bedding and has a significant impact on seismic wave propagation. Ignoring anisotropic effects during seismic migration degrades image quality. Therefore, we derive a pure qP-wave equation with high accuracy for modeling seismic wave propagation in tilted transversely isotropic (TTI) media. However, the derived pure qP-wave equation requires a computationally expensive spectral-based method for performing numerical simulations. This is unsuitable for large-scale industrial applications, particularly 3D applications. For numerical efficiency, we first decompose the newly derived wave equation into some fractional differential operators and spatial derivatives. The spatial derivatives can be directly solved by conventional finite-difference (FD) approaches. Then, we use an asymptotic approximation to find an equivalent form of fractional differential operators, obtaining scalar operators that we can discretize with the FD method. Numerical tests show that our TTI pure qP-wave equation with an FD discretization can produce accurate and highly efficient wavefield simulations in TTI media. We also use our TTI pure qP-wave equation with an FD discretization as a forward engine to implement TTI reverse time migration (RTM). Synthetic examples and a field data test demonstrate that our TTI RTM can effectively correct the anisotropic effects, providing high-quality imaging results while maintaining good computational efficiency.

Passive Source Reverse Time Migration Based on the Spectral Element Method

AbstractIncreasing deployment of dense arrays has facilitated detailed structure imaging for tectonic investigation, hazard assessment and resource exploration. Strong velocity heterogeneity and topographic changes have to be considered during passive source imaging. However, it is quite challenging for ray‐based methods, such as Kirchhoff migration or the widely used teleseismic receiver function, to handle these problems. In this study, we propose a 3‐D passive source reverse time migration strategy based on the spectral element method. It is realized by decomposing the time reversal full elastic wavefield into amplitude‐preserved vector P and S wavefields by solving the corresponding weak‐form solutions, followed by a dot‐product imaging condition to get images for the subsurface structures. It enables us to use regional 3‐D migration velocity models and take topographic variations into account, helping us to locate reflectors at more accurate positions than traditional 1‐D model‐based methods, like teleseismic receiver functions. Two synthetic tests are used to demonstrate the advantages of the proposed method to handle topographic variations and complex velocity heterogeneities. Furthermore, applications to the Laramie array data using both teleseismic P and S waves enable us to identify several south‐dipping structures beneath the Laramie basin in southeast Wyoming, which are interpreted as the Cheyenne Belt suture zone and agree with, and improve upon previous geological interpretations.

A Deep Learning Structure Correction Method Under the Influence of Special Lithology Body

Abstract In the process of precise exploration and development of oil and gas reservoirs, the accurate implementation of the target layer structure plays a vital role in the increase of oil reserves and the success of drilling. However, with the gradual deepening of the research, it is found that in the overlying strata of oil and gas reservoirs, special lithological bodies with strong heterogeneity and abnormal velocity are usually locally developed, which directly affects the construction of seismic data migration velocity, resulting in structural distortion at the local position of the reservoir in the process of data interpretation. Focusing on the structural change caused by the abnormal velocity, the deep learning algorithm is used to accurately predict the velocity of the special lithologic body and realize the accurate correction of the structure. Firstly, based on the combination of seismic and geological understanding, the geological model of the target layer and the overlying special lithologic body is constructed, and the distribution law of the special lithologic body is identified by the forward simulation method. Then use the well-side seismic trace and P-wave velocity curve to train the velocity data through the deep learning method and apply the training results to the post-stack seismic data to obtain the accurate P-wave velocity body. Finally, the real formation velocity over the target layer is restored to achieve reservoir structure correction. Based on this method, the research experiment was carried out in the M block of the Tarim Basin, and the spatial distribution range of the special lithology body was carved. The data error between the migration velocity body and the abnormal high-velocity body was calculated, and the target layer structure in the study area was corrected.

A microlocal and visual comparison of 2D Kirchhoff migration formulas in seismic imaging*

Abstract The term Kirchhoff migration refers to a collection of approximate linearized inversion formulas for solving the inverse problem of seismic tomography which entails reconstructing the Earth’s subsurface from reflected wave fields. A number of such formulas exists, the first dating from the 1950 s. As far as we know, these formulas have not yet been mathematically compared with respect to their imaging properties. This shortcoming is to be alleviated by the present work: we systematically discuss the advantages and disadvantages of the formulas in 2D from a microlocal point of view. To this end we consider the corresponding imaging operators in an unified framework as pseudodifferential or Fourier integral operators. Numerical examples illustrate the theoretical insights and allow a visual comparison of the different formulas.

Matrix imaging as a tool for high-resolution monitoring of deep volcanic plumbing systems with seismic noise

Volcanic eruptions necessitate precise monitoring of magma pressure and inflation for improved forecasting. Understanding deep magma storage is crucial for hazard assessment, yet imaging these systems is challenging due to complex heterogeneities that disrupt standard seismic migration techniques. Here we map the magmatic and hydrothermal system of the La Soufrière volcano in Guadeloupe by analyzing seismic noise data from a sparse geophone array under a matrix formalism. Seismic noise interferometry provides a reflection matrix containing the signature of echoes from deep heterogeneities. Using wave correlations resistant to disorder, matrix imaging successfully unscrambles wave distortions, revealing La Soufrière’s internal structure down to 10 km with 100 m resolution. This method surpasses the diffraction limit imposed by geophone array aperture, providing crucial data for modeling and high-resolution monitoring. We see matrix imaging as a revolutionary tool for understanding volcanic systems and enhancing observatories’ abilities to monitor dynamics and forecast eruptions.

Investigation of the 2011 Yingjiang, Yunnan, China Ms. 5.8 Earthquake Sequence: Seismic Migration, Seismogenic Mechanism, and Hazard Implication

On March 10, 2011, an Ms. 5.8 earthquake struck Yingjiang City, western Yunnan, China, causing destructive damage. Due to the very sparse distribution of seismic stations on the southwestern border of China, its seismogenic structure and mechanism remain controversial. In this study, with the aid of machine-learning-based detection and location workflow and template matching technique, we detect 10,356 events ranging from December 1, 2010, to April 30, 2011. The high-precision earthquake catalog shows that the foreshocks initiated in the extensional stepover connecting the northeast and middle segments of the Dayingjiang fault and then bilaterally extended northeast and southwest, with migration fronts that can be simulated by fluid diffusion model with diffusivities of 0.8 m2/s and 0.19 m2/s, respectively. The mainshock occurred at the southwest end of the foreshock sequence and then probably activated the northwest-trending blind fault. In addition, we determine the full moment tensor solutions for the mainshock, six large foreshocks, and one aftershock, with magnitudes ranging from 3.03 to 5.80, in which the mainshock was characterized by an obvious negative isotropic (ISO) component. The static Coulomb failure stress change caused by five Mw ≥ 4.0 foreshocks on the mainshock fault plane is ∼24 kPa, reaching the typical static triggering threshold. Therefore, we suggest that both the fluid diffusion and stress perturbation contribute to triggering the mainshock. This study advances our understanding of the spatiotemporal evolution, seismogenic mechanism, and hazard implication for the Yingjiang Ms. 5.8 earthquake and provides additional evidence of natural fluid-triggered seismicity in western Yunnan.

Seismic performance investigation of PCS water tank designed as tuned mass damper (TMD) for nuclear containment plant considering soil-structure interaction

Tuned mass damper (TMD) has been applied in civil engineering fields to reduce the seismic response. However, few TMDs are used in the nuclear structure due to the huge additional mass of the traditional TMD device on the nuclear structure. In this paper, the PCS water tank of a nuclear structure is reconstructed into a TMD system to reduce the seismic response of the nuclear structure. The proposed novel water tank-TMD device can use its characteristics of large mass to achieve the tuned shock absorption effectively while maintaining the original function of the water tank. Moreover, more and more nuclear power structures are being built on non-bedrock sites, and the effects of soil-structure interaction (SSI) on the seismic response and seismic migration performance of nuclear power structures need to be reasonably considered. The effects of SSI and auxiliary plants on the seismic response of the nuclear containment plants equipped with PCS water tanks are first investigated. Subsequently, the performance of the water tank-TMD system in reducing the seismic responses of nuclear containment plants is analyzed in detail by considering the effect of SSI, design frequency of TMD, and design criteria of TMD. The research results show that the proposed water tank-TMD system can effectively reduce the seismic response of nuclear containment plants. Besides, it is necessary to consider the influence of SSI on the dynamic characteristics of the structure in the design of TMD when the SSI effect is strong. The water tank-TMD system designed according to the TMD criteria proposed by Sadek et al., and Salvi and Rizzi has better reduction performance compared with that designed by other design criteria.

Angle-dependent image-domain least-squares migration through analytical point spread functions

Image-domain least-squares migration (IDLSM) is an established approach to recover high-fidelity seismic images of subsurface reflectors; this is achieved by removing the blurring effects of the Hessian operator in the standard migration approach with the help of so-called point spread functions (PSFs). However, most of the existing IDLSM approaches recover an angle-independent image of the subsurface reflectors, which is not suitable for subsequent amplitude-variation-with-angle (AVA) analysis. To overcome this limitation, we have developed an angle-dependent IDLSM approach, denoted as AD-IDLSM, which can recover a high-fidelity and high-resolution angle-dependent reflectivity image of subsurface reflectors. The problem is formulated here as an angle-dependent image-domain inversion with PSFs computed by means of full-wave Green’s function. More specifically, we derive an analytical expression to compute angle-dependent PSFs by means of a wave-equation-based Kirchhoff migration (WEBKM) engine, where a localization assumption is made in both spatial directions to decrease the computational cost and memory overhead. The amplitude and traveltime of the Green’s functions involved in the WEBKM approach are estimated by the excitation amplitude and excitation time of the full-wave wavefield. The scattering angle is then approximately estimated from the Poynting vector of the excitation-time field. To stabilize the solution of AD-IDLSM, we use a regularization scheme that applies a second derivative along the direction of the reflection angle of angle-domain common-image gathers (ADCIGs) to ensure continuity in the amplitude variations versus angle and suppress migration artifacts. We demonstrate the effectiveness of the AD-IDLSM approach through two synthetic and one field marine data set; the presented results confirm that AD-IDLSM can create ADCIGs with higher spatial resolution, better amplitude fidelity, and fewer migration artifacts compared with those obtained by its migration counterpart. Moreover, AD-IDLSM amplitude variations with angle are shown to closely resemble the theoretical AVA curve of the reflectors.

Neotectonics of the Barents Sea shelf eastern part: seismicity, faults and impact of the Atlantic-Arctic Rift System

distribution with the fault network defined from seismic data, and kinematic characteristics of seismic activity spatial migration are obtained. It is shown that seismic events recorded by the Norwegian NORSAR regional network within the Russian part of the Barents Sea shelf are grouped into linear clusters along shear kinematics faults. The fault network displaces Mesozoic seismic complexes and reaches the bottom surface, displacing quaternary deposits, which clearly indicates the modern age of displacements along which linear clusters of weak seismicity are grouped. The calculation of the total seismic moment in the space-time dimension showed the presence of seismic activity migration along short fragments of faults on the shelf with velocities from 10.5 to 25.7 km/year. There has been a surge in general activity in the shelf area since 2012. A comparison of the temporal evolution of seismic activity on the shelf with fragments of the Atlantic-Arctic Rift system suggests that it is the effect of tectonic deformation waves triggered along the geodynamically active intraplate boundary and propagating to the shelf at a speed of 20–22 km/year. Migration rates with speeds up to 77 km/year are less likely. There is a possibility that the increase in the intensity of seismic activity on the shelf after 2012 is not an emission from the effects of a slow deformation wave, but the result of a direct trigger effect on the shelf from the structures of the Knipovich and Gakkel ridges.

3D reverse-time migration for pure P-wave in orthorhombic media

Compared with the transverse isotropic (TI) medium, the orthorhombic anisotropic medium has both horizontal and vertical symmetry axes and it can be approximated as a set of vertical fissures developed in a group of horizontal strata. Although the full-elastic orthorhombic anisotropic wave equation can accurately simulate seismic wave propagation in the underground media, a huge computational cost is required in seismic modeling, migration, and inversion. The conventional coupled pseudo-acoustic wave equations based on acoustic approximation can be used to significantly reduce the cost of calculation. However, these equations usually suffer from unwanted shear wave artifacts during wave propagation, and the presence of these artifacts can significantly degrade the imaging quality. To solve these problems, we derived a new pure P-wave equation for orthorhombic media that eliminates shear wave artifacts while compromising computational efficiency and accuracy. In addition, the derived equation involves pseudo-differential operators and it must be solved by 3D FFT algorithms. In order to reduce the number of 3D FFT, we utilized the finite difference and pseudo-spectral methods to conduct 3D forward modeling. Furthermore, we simplified the equation by using elliptic approximation and implemented 3D reverse-time migration (RTM). Forward modeling tests on several homogeneous and heterogeneous models confirm that the accuracy of the new equation is better than that of conventional methods. 3D RTM imaging tests on three-layer and SEG/EAGE 3D salt models confirm that the ORT media have better imaging quality.

Seismic migration of water‐bottom‐related multiples accelerated by random phase‐encoding strategy

AbstractMarine seismic multiples contain more structure information than primaries and should be considered in migrations. However, multiple migrations suffer from severe crosstalks generated by interferences among undesirable multiples. It has been proven that the water‐bottom‐related multiple migration can suppress crosstalks greatly. However, if all associated consecutive‐order multiples are considered, the computation cost is extremely high. To settle this issue, a phase‐encoding‐based multiple migration is proposed. Supergathers are first created by randomly phase‐encoding consecutive‐order multiples and stacking‐encoded multiples. By migrating supergathers, the proposed method can fulfil migrations of all order multiples simultaneously, thereby reducing the computation cost significantly. We use a three‐layer model and the Pluto 1.5 model for numerical comparisons. The results reveal that the method can retrieve high‐quality images and increase computation efficiency considerably.

High-precision imaging of small voids in tunnel lining based on reverse time migration

Abstract Diseases such as voids behind the initial support of the tunnel will lead to problems such as lining rupture and concrete damage in the tunnel structure, which seriously affects the driving safety in the tunnel. In tunnel construction, ground penetrating radar is frequently employed as a method for detecting hidden defects. However, when the cavity size is small, there will be a large error when the original radar data is directly interpreted, which cannot meet the needs of practical engineering. To enhance the precision in tunnel detection, the Finite-Difference Time-Domain method is used to simulate various small cavities of different shapes located behind the primary support lining of the tunnel, and the study involves examining the electromagnetic response characteristics of different kinds of holes. The migration techniques of KIRCHHOFF, F-K, and reverse time are employed to reconstruct the hole target. Furthermore, digital morphology is employed to enhance the clarity of the image. Finally, the edge detection technology is used to further extract the hole features. After that, small cavities of different shapes are buried in the established outdoor concrete model box and radar detection is carried out. The detection position of the reconstructed radar image and the actual position of the cavity in the model box are verified. The results show that compared with KIRCHHOFF migration and F-K migration, the image processed by reverse time migration is clearer and more intuitive, and can restore the contour of small holes well. The research findings can serve as a reference for the interpretation of radar data of the cavity behind the initial support of the tunnel.

Depth imaging of multicomponent seismic data through the application of 2D full‐waveform inversion to P‐ and SH‐wave data: SEAM II Barrett model study

AbstractWe examine the value of the nine‐component seismic survey by generating the Kirchhoff depth migration images of compressional wave (P‐wave), converted wave (PS‐wave) and horizontally polarized shear wave (SH shear wave) data simulated from the SEAM II Barrett unconventional model. We first utilize full waveform inversion to obtain a P‐wave velocity model from P‐wave data and an S‐wave velocity model from SH‐wave data. Both P‐wave and SH‐wave data are generated with the maximum frequency of 20 Hz while assuming that the subsurface is isotropic. To implement full waveform inversion, we use a two‐dimensional time‐domain finite‐difference method and the L2 norm to measure the data misfit. We use both refractions and reflections in P‐ and SH‐wave data to reconstruct the P‐ and S‐wave velocity models from the surface to the reservoir. The inverted P‐ and S‐wave velocities contain the main features of the model (e.g., major faults and channels) but have some difficulties in estimating high‐frequency velocity variation within the first 300‐m depth of the model due to the frequency constraint. We then use the inverted P‐ and S‐wave velocities to generate Kirchhoff depth migration gathers and images from the P‐, PS‐ and SH‐wave data. The flat P‐ and SH‐wave common‐image offset gathers suggest that SH‐ and P‐wave full waveform inversion can generate adequate S‐ and P‐wave velocities for migration. Flat PS‐wave gathers and the clear PS‐wave migration image are also obtained using the inverted P‐ and S‐wave velocities simultaneously. This result indicates that obtaining S‐wave velocities from SH‐wave data can aid PS‐wave data processing and imaging. Moreover, the SH‐wave images and S‐wave images of the radial component provide better delineation of fault planes and small‐scale geobodies within the reservoir since the wavelength of the S‐wave is smaller compared to P‐wave when similar frequency ranges are recorded. Therefore, our study shows that S‐wave velocities can be successfully constructed by the two‐dimensional full waveform inversion application of the SH‐wave data. The subsequent imaging of multicomponent seismic data improves the delineation of certain unconventional reservoirs compared to the traditional P‐wave imaging.

A shortest‐path‐aided fast‐sweeping method to improve the accuracy of traveltime calculation in vertically transverse isotropic media

AbstractThe high accuracy and efficiency of traveltime calculation are critical in seismic tomography, migration, static corrections, source locations and anisotropic parameter estimation. The fast‐sweeping method is an efficient upwind finite‐difference approach for solving the eikonal equation. However, the fast‐sweeping method is accurate only along the axis directions. In two‐dimensional or higher dimensional cases, the accuracy is severely decreased in the diagonal directions due to the numerical errors in these directions. These similar numerical errors also arose in higher order fast‐sweeping method and anisotropic fast‐sweeping method. To improve the accuracy of traveltime calculation in two‐dimensional or higher dimensional space, a shortest‐path‐aided fast‐sweeping method is proposed. The shortest‐path‐aided solution is embedded into the sweeping process of the standard fast‐sweeping method to improve the traveltime accuracy in the diagonal directions. Shortest‐path‐aided fast‐sweeping method is very easy to implement nearly without additional computational cost and memory consumption. Furthermore, this method is easy to extend from two‐dimensional to higher dimensional, from low‐order to higher‐order and from isotropic to anisotropic cases.

Distributed acoustic sensing microseismic reflection imaging for hydraulic fracture and fault lineament characterization

This study presents a novel workflow designed for migrating reflected S waves generated by microseismic events, as recorded by downhole distributed acoustic sensing (DAS), to characterize hydraulic fractures in three dimensions. In contrast to existing fracture imaging techniques, which have encountered challenges in accurately representing fracture networks and often rely on simplified models, our imaging technique does not assume that fractures are planar or in a prespecified orientation. DAS seismic measurements benefit from the large aperture and dense spatial sampling enabled by the kilometers-long fiber and, therefore, are able to capture a large number of strong reflections compared with traditional borehole geophones or accelerometers. We treat microseismic events as high-frequency sources and apply prestack Kirchhoff migration to each individual source after wavefield separation. Fracture imaging results for multiple selected events are then stacked to generate a 3D reflectivity volume, revealing the subsurface fracture and fault networks in intricate detail. The high-resolution fracture images generated by the developed reflection migrating process illuminate the heart of the stimulated volume of the reservoir, a zone that is often challenging to access using conventional surface arrays or active sources. To validate the effectiveness of our workflow, our study uses a data set acquired during a multiwell project in the Eagle Ford Shale and Austin Chalk in South Texas. To assess the accuracy and reliability of the results, the reflection imaging output is integrated with the microseismicity distribution and strain measurements from low-frequency DAS for interpretation. The results of reflection imaging improve our understanding of fracture geometry, including distal fractures that are away from the monitoring well, allow the direct estimation of fracture height and length, and potentially signify the presence of preexisting fluid-filled fault lineaments.

Adaptive structure-based full-waveform inversion

In simultaneous-shot full-waveform inversion (FWI), the re-started limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm can be used to suppress crosstalk, but crosstalk cannot be completely eliminated from the inversion results. To further solve the crosstalk problem caused by the interference among individual shots in simultaneous-shot FWI, adaptive structure-based smoothing is applied to the FWI with the restarted L-BFGS algorithm. Structure-based smoothing can mitigate the crosstalk by highlighting structure boundaries. To perform structure-based smoothing, an implicit anisotropic diffusion equation is solved. We carry out a multiscale inversion strategy. Because the FWI results in the low-frequency band contain less structural information and more crosstalk, structure-based smoothing is applied to the frequency band at frequencies higher than the peak frequency to prevent negative effects from the model structure in the low-frequency band. The estimated structural information is iteratively updated during the inversion process. Furthermore, structure-based smoothing is only added to the iterations with invariant encoding to reduce oversmoothing. The invariant encoding means that the shot encoding remains unchanged between iterations. Numerical experiments with an overthrust model indicate that the proposed adaptive structure-based FWI method can effectively suppress artifacts and provide a high-resolution inversion result, even when the encoded data are contaminated with strong noise. Another advantage of the adaptive structure-based FWI method is that no seismic migration needs to be performed, which makes it more efficient for real data.

Analysis of Seismic Response Characteristics of Fractured Carbonate Reservoirs Based on Physical Model (Tarim Basin)

Anomalous bright spots, called the string of bead-like response, are typical seismic migration profile features in carbonate fractured reservoirs in the Tarim Basin, and they are indicators of high-quality oil and gas reservoirs. Correctly recognizing the correspondence between fractures and the SBLR can contribute to the efficient drilling of target carbonate fractured reservoirs. Physical models can describe fractured reservoirs more directly and efficiently than real situations and have obvious advantages in accurately and quantitatively designing parameters such as dipping angles and the number of layers of fractured reservoirs. Under such a background, according to the real tectonic characteristics of the Tarim Basin, among the main hydrocarbon reservoirs, fractured reservoirs with various parameters were designed and a physical model was constructed according to the real stratigraphic parameters. After seismic data acquisition and processing, the response characteristics of the string of bead-like response were extracted and summarized from seismic migration profiles for all fractured reservoirs, which provided targeted analyses and discussions on the fracture dipping angle, number of fracture layers, overlying stratigraphic influences, and planar attributes of the fractured reservoirs. In general, the larger the fracture dip, the more difficult it is to identify, while the slope of reflection strength and maximum absolute amplitude attributes can be important markers for fractured reservoir identification. The physical modeling study of fractured reservoirs in this paper can provide a basis for the analysis and prediction of carbonate fractured reservoirs in the Tarim Basin.

DiffraPy: An open-source Python software for seismic diffraction imaging

Seismic diffraction events, often overlooked as noise in conventional seismic imaging methods, contain valuable high-resolution information about subsurface structures and heterogeneities. Furthermore, the diffractive component of the wavefield is masked by the stronger reflection events. Although the literature presents innovative methods for separating diffracted events from the total wavefield, there are currently limited options of open-source alternatives in seismic diffraction processing and imaging. In this study, we introduce DiffraPy, an open-source Python package specifically designed for performing both conventional and diffraction-oriented migrations. The main workflow involves implementing an anti-stationary phase filter in the Kirchhoff migration for enhancing energy outside the Fresnel zone. Additionally, we propose using the semblance matrix obtained from the dip field calculation as an additional weight in conventional migration. The code provides functionality for forward modeling seismic data, allowing for new testing and experimentation. We demonstrate the general workflow of the code using a toy model and evaluate its capabilities on a synthetic salt model. In both cases, the program successfully suppressed the majority of reflected energy in the diffraction-oriented image. The final image highlights the responses of small lenses, discontinuities, and faults, providing valuable insights into the subsurface for interpretation. Despite the method sensitivity to variations in the correct velocity model, we demonstrate that our code is capable of imaging diffractors even when utilizing incorrect migration velocities of up to 5%. Our wavenumber analysis reveals the diffraction image preservation of smaller wavelengths, indicating higher resolution compared to conventional products. Future implementations may include extension to 3D datasets and improve the computational efficiency of our code. We expect our code to offer the geosciences community the first freely accessible and well-documented Python alternative for testing, reproducing, and exploring diffraction imaging examples.