Wave Equation Migration Methods Research Articles

Reverse time migration surface offset gathers by attribute migration and constrained least-squares inversion

Abstract Surface offset gathers (SOGs) are crucial for velocity updating in seismic data migration and muting the stretched waveform in shallow parts of migrated seismic data. Kirchhoff-based migration methods can output SOGs because they take every seismic trace as the input. As a result, the imaging results can be rearranged according to the offset value of the input data to obtain SOGs. However, regarding more accurate wave equation methods (such as reverse time migration, RTM), the SOGs cannot be obtained directly with a single migration calculation. Attribute migration is an easy-to-implement method to output SOGs of wave equation migration methods. It calculates the offset value of each imaging point by dividing the migration result modulated by the offset attribute with the conventional migration result. However, the division error may cause the calculation result outside the given offset value. This paper proposes a constrained least-squares inversion method for stable calculation of the RTM offset range. The method ensures that the calculated results are within the given offset range. We tested the method against the direct division and least-squares methods without constraints using the Marmousi model and a real 3D dataset. We show that the SOGs obtained by the constrained least-squares inversion attribute migration method had the same characteristics as those derived using the division-based attribute migration method, but they also had higher vertical resolution, better energy distribution, and improved lateral continuity. In a further study, we expect to develop a stable and direct method to efficiently calculate SOGs for RTM and an iterative closed loop for RTM velocity updating.

Lithospheric structure beneath the boundary region of North China Craton and Xing Meng Orogenic Belt from S-receiver function analysis

The boundary region of the North China Craton (NCC) and Xing-Meng Orogenic Belt (XMOB), including the northern boundary of the NCC and the southern XMOB, is the ideal area to investigate the spatially uneven lithospheric deformation and associated tectonic evolution in the boundary region of the craton and orogen. The depths and structures of the lithosphere–asthenosphere boundary (LAB) and mid-lithosphere discontinuity (MLD) are the key information to describe the nature and deformation pattern of the lithosphere, and thus can provide basic constraints on the associated tectonic processes and mechanism of the lithospheric deformation. Based on waveform data collected from 226 broadband seismic stations, high resolution spatial variations in the LAB and MLD depths were obtained by using the wave equation migration method for the S-receiver functions. The migrated images show a lateral variation in the LAB depth from 100–120 km in the northern boundary of the NCC to 120–140 km in the southern XMOB. The imaged LAB within the NCC is much shallower than that revealed by mantle xenoliths enclosed in the nearby Early Paleozoic kimberlites (~180 km), suggesting that the cratonic lithosphere of the area may have undergone significant thinning during the Phanerozoic evolution. Furthermore, a coherent MLD is identified at depths of 70–90 km beneath the northeastern NCC, which corroborates the speculation that the MLD probably existed in the eastern NCC before lithospheric destruction in the Mesozoic. The spatial variations in the LAB and MLD depths probably reflect complex lithospheric deformation and thinning around the boundary region of NCC and XMOB. The properties, including the gradient thickness and S-wave velocity contrast, of the two discontinuities further indicate that lithospheric deformation and thinning might be related to partial melting induced by the Big Mantle Wedge associated with the subducting Paleo-Pacific plate and its stagnation during the Late Mesozoic.

Amplitude-Preserved Wave Equation: An Example to Image the Gas Hydrate System

Natural gas hydrate is an important energy source. Therefore, it is extremely important to provide a clear imaging profile to determine its distribution for energy exploration. In view of the problems existing in conventional migration methods, e.g., the limited imaging angles, we proposed to utilize an amplitude-preserved one-way wave equation migration based on matrix decomposition to deal with primary and multiple waves. With respect to seismic data gathered at the Chilean continental margin, a conventional processing flow to obtain seismic records with a high signal-to-noise ratio is introduced. Then, the imaging results of the conventional and amplitude-preserved one-way wave equation migration methods based on primary waves are compared, to demonstrate the necessity of implementing amplitude-preserving migration. Moreover, a simple two-layer model is imaged by using primary and multiple waves, which proves the superiority of multiple waves in imaging compared with primary waves and lays the foundation for further application. For the real data, the imaging sections of primary and multiple waves are compared. We found that multiple waves are able to provide a wider imaging illumination while primary waves fail to illuminate, especially for the imaging of bottom simulating reflections (BSRs), because multiple waves have a longer travelling path and carry more information. By imaging the actual seismic data, we can make a conclusion that the imaging result generated by multiple waves can be viewed as a supplementary for the imaging result of primary waves, and it has some guiding values for further hydrate and in general shallow gas exploration.

Least-squares Kirchhoff migration with non-smooth regularisation strategy for subsurface imaging

During the past several decades, many types of wave-equation migration methods arise for subsurface structure imaging. The classical Kirchhoff migration, however, is still widely adopted in the petroleum industry owing to its flexibility and computational efficiency. In constant density isotropic acoustic media, a basic assumption of the Kirchhoff migration is that every point of the subsurface model is supposed to be a diffractor which scatters wavefield energy to every direction, and hence collecting the scattered energy of all directions is the basic requirement for focusing the diffractor. Factors influencing the final image quality include incomplete data acquisition, multipathing from the surface to the imaging point, and insufficient illumination under complex overburden. All these factors can be theoretically taken into account in the migration weighting coefficient. However, computation of the weighting coefficient is hard work. In view of this difficulty, a fast regularising least-squares Kirchhoff migration algorithm is presented in this paper. It not only accounts for the irregular and incomplete data sampling (e.g. limited recording aperture, coarse sampling and acquisition gaps), but also compensates for the anomalous ray coverage and multipathing problem except for the shadow zone in the media. For the purpose of attenuating migration artefacts and providing a clear and accurate image of subsurface reflectivity, regularisation strategies are applied. The classical regularisation strategy may easily lead to over-regularisation or insufficient regularisation; we try to balance these two effects in this paper. The method is called the hybrid regularisation which incorporates smoothing and non-smoothing scale operators. The algorithm is implemented using a fast gradient decent solution method based on the Rayleigh quotient being used. Numerical experiments show that this hybrid regularisation method is powerful in handling the sparsity and smoothness of the model parameters.

Superresolution imaging using resonant multiples

A resonant multiple is defined as a multiple reflection that revisits the same subsurface location along coincident reflection raypaths. We have found that resonant first-order multiples can be migrated with either Kirchhoff or wave-equation migration methods to give images with approximately twice the spatial resolution compared with poststack primary-reflection images. We have developed a moveout-correction stacking method to enhance the signal-to-noise ratios of the resonant multiples for superresolution migration. The effectiveness of this procedure is validated by synthetic and field data tests.

Importance of Seismic Diffractions for Fractures Imaging

The use of diffraction imaging is important and rapidly gaining momentum in the oil and gas industry as the need of the industry moves toward exploiting smaller and more complex structures to find hydrocarbon. Illumination of these small scale discrete transmissivity structures such as faults or fracture zones prior to exploration offers substantial benefits for all phases of field development. We have established a new way of seismic imaging through diffractions study using reflection seismology. A subsurface geological model is developed from one of highly productive field located in the southern part of the Malay Basin. To study the seismic component in reflection seismology, a Finite difference modelling is used to generate synthetic seismic data. Both velocity and density models are used for the explanation of wave propagation, intimating the subsurface from side to side in both one-way and two-wave wave propagation. One–way wave equation migration methods are tested to image the faults in the synthetic seismic section of the Malay basin and we demonstrated the need of diffraction imaging for highly dipping faults and complex structure.

Kinematic artifacts in the subsurface-offset extended image and their elimination by a dip-domain specularity filter

Common-image gathers in the dip-angle domain may be computed in relation to wave-equation migration methods, extended by the subsurface offset. They involve the application of a postmigration local Radon transform on the subsurface-offset extended image. In the dip-angle domain, seismic reflections are focused around the specular dip angle of reflection. This focusing distinguishes them from any other event in the image space. We have incorporated the dip-angle information about the presence of specular reflections into the computation of the conventional scattering-angle-dependent reflection coefficient. We have designed a specularity filter in the dip-angle domain based on a local semblance formula that recognizes and passes events associated with specular reflections, while suppressing other sorts of nonspecular signal. The filter is remarkably effective at eliminating either random or coherent noises that contaminates the prestack image. In particular, our dip-angle filter provides a method for the suppression of kinematic artifacts, commonly generated by migration in the subsurface-offset domain. These artifacts are due to an abrupt truncation of the data acquisition geometry on the recording surface. We have studied their appearance and devised an appropriate formation mechanism in the subsurface-offset and scattering-angle domains. The prominent presence of the kinematic artifacts in image gathers usually impairs the quality of the postmigration analysis and decelerates the convergence of wave-equation inversion techniques. We have determined from testing on synthetic and field data that using the proposed dip-angle-domain specularity filter efficiently eliminates the kinematic artifacts in the delivered gathers. We expect involvement of the specularity filter to increase the reliability and quality of the seismic processing chain and provide a faster convergence of iterative methods for seismic inversion.

Scattering and dip angle decomposition based on subsurface offset extended wave-equation migration

An angle-dependent reflection coefficient is recovered by seismic migration in the angle domain. We have developed a postmigration technique for computing scattering and dip angle common-image gathers (CIGs) from seismic images, extended by the subsurface offset, based on wave-equation migration methods. Our methodology suggests a system of Radon transform operators by introducing local transform relations between the subsurface offset image and the angle-domain components. In addition to the commonly used decomposition of the scattering angle, the methodology associates the wave-equation migration with dip-domain images as well. The same postmigration subsurface offset image is used to decompose scattering and dip angle CIGs individually or to decompose a multiangle CIG by showing simultaneously both angles on the gather’s axis. We show that the dip-angle response of seismic reflections is a spot-like signature, focused at the specular dip of the subsurface reflector. It differs from the well-studied smile-like response usually associated with reflections in the dip domain. The contradiction is clarified by the nature of the subsurface offset extension, and by emphasizing that the angles are decomposed from the subsurface offset image after the imaging condition, without directly involving the propagating incident and scattered wavefields. Several synthetic and field data tests proved the robustness of our decomposition technique, by handling various subsurface models, including seismic diffractions. It is our belief that dip-angle information, decomposed by wave-equation migration, would have a great impact in making the scattering-angle reflection coefficient more reliable and noise free, in addition to the acceleration of wave-equation inversion methods.

Near-Surface Multi-Point Vibration Location Method Based on Seismic Depth Migration

This paper presents a method for multi-point vibration location near-surface. Unlike traditional source location technologies, it is not using travel times of seismic waves for positioning calculation, but according to the entire seismic data record, using the wave equation migration method to calculate the source location. Similar to exploration seismic, the records from a survey line within a certain period of time are data volumes with dimensions of time and ground coordinates. Based on the data, combined with surface seismic wave propagation characteristics, by using the improved seismic depth migration algorithm, the vibration energy will return to the starting position where exactly the source location is. The method can solve the problem of location calculation error by using traditional method when several vibration at the same time or continuous vibration occurs at some point.

An effective denoising strategy for wave equation migration based on propagation angles

We present an effective denoising strategy for two-way wave equation migration. Three dominant artifact types are analyzed and eliminated by an optimized imaging condition. We discuss a previously unsolved beam-like artifact, which is probably caused by the cross-correlation of downward transmitting and upward scattering waves from both the source and receiver side of a single seismic shot. This artifact has relatively strong crosscorrelation but carries no useful information from reflectors. The beam-like artifact widely exists in pre-stack imaging and has approximately the same amplitude as useful seismic signals. In most cases, coherent artifacts in the image are caused by directionally propagating energy. Based on propagation angles obtained by wavefield gradients, we identify the artifact energy and subtract its contribution in the imaging condition. By this process most artifacts can be accurately eliminated, including direct wave artifacts, scattering artifacts, and beamlike artifacts. This method is independent of the wavefield propagator and is easy to adapt to almost all current wave equation migration methods if needed. As this method deals with the physical artifact origins, little damage is caused to the seismic signal. Extra k-domain filtering can additionally enhance the stacking result image quality. This method succeeds in the super-wide-angle one-way migration and we can expect its success in other two-way wave equation migrations and especially in reverse time migration.

3D angle gathers from reverse time migration

Common-image gathers are an important output of prestack depth migration. They provide information needed for velocity model building and amplitude and phase information for subsurface attribute interpretation. Conventionally, common-image gathers are computed using Kirchhoff migration on common-offset/azimuth data volumes. When geologic structures are complex and strong contrasts exist in the velocity model, the complicated wave behaviors will create migration artifacts in the image gathers. As long as the gather output traces are indexed by any surface attribute, such as source location, receiver location, or surface plane-wave direction, they suffer from the migration artifacts caused by multiple raypaths. These problems have been addressed in a significant amount of work, resulting in common-image gathers computed in the reflection angle domain, whose traces are indexed by the subsurface reflection angle and/or the subsurface azimuth angle. Most of these efforts have concentrated on Kirchhoff and one-way wave-equation migration methods. For reverse time migration, subsurface angle gathers can be produced using the same approach as that used for one-way wave-equation migration. However, these approaches need to be revisited when producing high-quality subsurface angle gathers in three dimensions (reflection angle/azimuth angle), especially for wide-azimuth data. We have developed a method for obtaining 3D subsurface reflection angle/azimuth angle common-image gathers specifically for the amplitude-preserved reverse time migration. The method builds image gathers with a high-dimensional convolution of wavefields in the wavenumber domain. We have found a windowed antileakage Fourier transform method that leads to an efficient and practical implementation. This approach has generated high-resolution angle-domain gathers on synthetic 2.5D data and 3D wide-azimuth real data.

The Algorithm of High Order Finite Difference Pre‐Stack Reverse Time Migration and GPU Implementation

AbstractPre‐stack reverse time migration (RTM) is a very useful tool for seismic imaging. But it has not been widely used because of the highly intensive computation cost, imaging noise and massy memory demand. In this paper, we present the implementation process of RTM and analyze the stability condition and dispersion relation of the finite difference (FD) method. For the problem of intensive computation cost, we use the Graphic Processing Unit (GPU) architecture to realize RTM and get an order of magnitude higher speedup ratio compared to the traditional CPU architecture. Compared to the one way wave equation migration methods, RTM does not have the imaging dip limit and the imaging effect is significantly improved. The tests on vertical fault and BP models prove this conclusion. The problems of imaging noise removal and massy memory demand will be discussed in other papers.

Overview and classification of wavefield seismic imaging methods

The literature and seismic processing practice overflow with a diversity of wave-equation migration methods. So what makes wave-equation migration a wave-equation migration? What is the commonality of these methods and how do they differ? This article provides a framework for understanding different wave-equation migration methods.

The finite-frequency sensitivity kernel for migration residual moveout and its applications in migration velocity analysis

We have derived a broadband sensitivity kernel that relates the residual moveout (RMO) in prestack depth migration (PSDM) to velocity perturbations in the migration-velocity model. We have compared the kernel with the RMO directly measured from the migration image. The consistency between the sensitivity kernel and the measured sensitivity map validates the theory and the numerical implementation. Based on this broadband sensitivity kernel, we propose a new tomography method for migration-velocity analysis and updating — specifically, for the shot-record PSDM and shot-index common-image gather. As a result, time-consuming angle-domain analysis is not required. We use a fast one-way propagator and multiple forward scattering and single backscattering approximations to calculate the sensitivity kernel. Using synthetic data sets, we can successfully invert velocity perturbations from the migration RMO. This wave-equation-based method naturally incorporates the wave phenomena and is best teamed with the wave-equation migration method for velocity analysis. In addition, the new method maintains the simplicity of the ray-based velocity analysis method, with the more accurate sensitivity kernels replacing the rays.

Inversion of Ground‐Penetrating Radar Data for 2D Electric Parameters

AbstractThis paper derives inversion formulas for reconstructing electric parameters of a medium from observed ground‐penetrating radar (GPR) data based on the two‐dimensional (2D) Maxwell's equations. By introducing an adjoint field calculated using the Maxwell's equations, the gradients of the objective function with respect to unknown parameters can be expressed explicitly. The procedure of inversion is, first of all, to build an initial model and to calculate guessed data. Secondly, to calculate adjoint fields by solving the adjoint wave equations, and thirdly to calculate the gradients and to modify the model. These steps are repeated until an acceptable result is obtained. In the modification of the model, the conjugate gradient method is used. Three reconstruction examples are illustrated. A permittivity profile is reconstructed from synthetic data received on a rough surface. A permittivity and a conductivity profile are reconstructed simultaneously from experimental data. Reconstruction results obtained by using the method of this paper are compared with a result obtained by using the inversion method based on the scalar wave equation where absorption is neglected, and with a migration profile obtained by using a wave equation migration method.

A wave equation migration method for receiver function imaging: 1. Theory

A wave equation‐based poststack depth migration method is proposed to image the Earth's internal structure using teleseismic receiver functions. By utilizing a frequency wave number domain one‐way phase screen propagator for wave field extrapolation in the migration scheme, common conversion point (CCP) stacked receiver functions are backward propagated to construct a subsurface structural image. The phase screen propagator migration method takes into account the effects of diffraction, scattering, and travel time alternation caused by lateral heterogeneities, and therefore it is particularly useful for imaging complex structures and deep discontinuities overlain by strong shallow anomalies. Synthetic experiments demonstrate the validity of the migration method for a variety of laterally heterogeneous models. The migrated images show considerable improvement over the CCP images in recovering the model features. Influences of several factors on the image quality of the poststack migration are further investigated, including interstation spacing, noise level of the data, velocity model used in migration, and earthquake distribution (incident direction of source fields). On the basis of the sampling theorem and previous statistic results, we discuss the relation of spatial resolution and signal‐to‐noise ratio of the migrated image with the frequency of the data, surface station spacing and number of receiver functions used in stacking. We show that both CCP stacking and poststack migration of receiver functions need to be designed in a target‐oriented way for reliable and efficient imaging. Our results also suggest that careful consideration of earthquake source distribution is necessary in designing seismic experiments aimed at imaging steeply dipping structures.

A wave equation migration method for receiver function imaging: 2. Application to the Japan subduction zone

The newly developed wave equation poststack depth migration method for receiver function imaging is applied to study the subsurface structures of the Japan subduction zone using the Fundamental Research on Earthquakes and Earth's Interior Anomalies (FREESIA) broadband data. Three profiles are chosen in the subsurface imaging, two in northeast (NE) Japan to study the subducting Pacific plate and one in southwest (SW) Japan to study the Philippine Sea plate. The descending Pacific plate in NE Japan is well imaged within a depth range of 50–150 km. The slab image exhibits a little more steeply dipping angle (∼32°) in the south than in the north (∼27°), although the general characteristics between the two profiles in NE Japan are similar. The imaged Philippine Sea plate in eastern SW Japan, in contrast, exhibits a much shallower subduction angle (∼19°) and is only identifiable at the uppermost depths of no more than 60 km. Synthetic tests indicate that the top 150 km of the migrated images of the Pacific plate is well resolved by our seismic data, but the resolution of deep part of the slab images becomes poor due to the limited data coverage. Synthetic tests also suggest that the breakdown of the Philippine Sea plate at shallow depths reflects the real structural features of the subduction zone, rather than caused by insufficient coverage of data. Comparative studies on both synthetics and real data images show the possibility of retrieval of fine‐scale structures from high‐frequency contributions if high‐frequency noise can be effectively suppressed and a small bin size can be used in future studies. The derived slab geometry and image feature also appear to have relatively weak dependence on overlying velocity structure. The observed seismicity in the region confirms the geometries inferred from the migrated images for both subducting plates. Moreover, the deep extent of the Pacific plate image and the shallow breakdown of the Philippine Sea plate image are observed to correlate well with the depth extent of the seismicity beneath NE and SW Japan. Such a correlation supports the inference that the specific appearance of slabs and intermediate‐depth earthquakes are a consequence of temperature‐dependent dehydration induced metamorphism occurring in the hydrated descending oceanic crust.

High-resolution wave-equation AVA imaging: Algorithm and tests with a data set from the Western Canadian Sedimentary Basin

This paper presents a 3D least-squares wave-equation migration method that yields regularized common-image gathers (CIGs) for amplitude-versus-angle (AVA) analysis. In least-squares migration, we pose seismic imaging as a linear inverse problem; this provides at least two advantages. First, we are able to incorporate model-space weighting operators that improve the amplitude fidelity of CIGs. Second, the influence of improperly sampled data (footprint noise) can be diminished by incorporating data-space weighting operators. To investigate the viability of this class of methods for oil and gas exploration, we test the algorithm with a real-data example from the Western Canadian Sedimentary Basin. To make our problem computationally feasible, we utilize the 3D common-azimuth approximation in the migration algorithm. The inversion algorithm uses the method of conjugate gradients with the addition of a ray-parameter-dependent smoothing constraint that minimizes sampling and aperture artifacts. We show that more robust AVA attributes can be obtained by properly selecting the model and data-space regularization operators. The algorithm is implemented in conjunction with a preconditioning strategy to accelerate convergence. Posing the migration problem as an inverse problem leads to enhanced event continuity in CIGs and, hence, more reliable AVA estimates. The vertical resolution of the inverted image also improves as a consequence of increased coherence in CIGs and, in addition, by implicitly introducing migration deconvolution in the inversion.

Seismic migration on the IBM 3090 Vector Facility

Seismic prospecting aims at determining the structure of the earth from indirect measurements. Acoustic wave fields are generated at the surface, penetrate the earth, and are backscattered by the earth's inhomogeneities. The data recorded at the surface are processed in a complex sequence of steps among which seismic migration plays an important role. This is a wave depropagation process that permits the localization in depth of the origin of the diffraction events measured (in time) at the surface. This paper presents an overview of the major wave-equation migration methods. The most frequently executed algorithms or kernels on which the execution speed depends most crucially are given particular attention. The speedup resulting from scalar-to-vector formulation is presented over wide ranges of dimensionality for linear tridiagonal equation solvers, Fourier Transforms, and convolution operations. The vectorizability and resulting speedup are also addressed in the case of migration schemes known as the Phase-Shift Method and the Phase Shift Plus Interpolation (PSPI) Method. It is shown that Fourier domain migration based on the phase-shift concept lends itself conveniently to multilevel parallelism on the 3090 Vector Facility (VF): vectorization of the innermost loops and concurrent processing in the outer loops by means of the VS FORTRAN Version 2 Multitasking Facility.