Wave-equation Migration Research Articles

One-step sparse-matrix method for full wave equation depth migration

By completely solving the second-order full-acoustic wave equation, two-way/full wave equation migration generates images up to arbitrary angles in inhomogeneous media. Full wave equation depth migration (FWDM) necessitates less wavefield storage during imaging, generates clearer imaging results, and has more computational flexibility compared with conventional reverse time migration. Conventionally, FWDM uses Runge-Kutta or two explicit one-way operators to extrapolate the data, which involves either a multirecursive process or complex square root computations of the Helmholtz operator, causing implementation complexities and inconveniences. To overcome these challenges, a novel algorithm called the one-step sparse-matrix method for FWDM (OSM-FWDM) is developed and executed. This method involves a straightforward calculation of the full wave equation propagator through a Taylor expansion and using the so-called precise integration method. The computational efficiency can be significantly improved by using sparse-matrix multiplication in the implementation of wavefield extrapolation. Building upon that, the method serves as a basic solution complemented by FWDM, offering straightforward, efficient, and easy coding characteristics. Our OSM-FWDM method is simpler and faster compared with conventional FWDM, all while maintaining imaging accuracy. In addition, a new solution that uses an extended velocity model to remove source-related artifacts in FWDM is developed and illustrated. These improvements are demonstrated through a series of numerical examples and a real data application.

Image-domain least-squares imaging in an elastic vertically transversely isotropic medium: Multiparameter point-spread function deconvolution

By solving a linear inverse problem, least-squares migration (LSM) can provide higher-precision images of complex subsurface structures than traditional adjoint-operator-based migration. LSM can be conducted in the data domain and the image domain (IDLSM). IDLSM in acoustic media has been extensively studied, commonly using single-parameter point-spread function (PSF) deconvolution. For high-accuracy imaging of complex media, we develop an image-domain elastic least-squares imaging method using multiparameter PSF deconvolution for vertically transversely isotropic (VTI) media. The multiparameter PSFs, including P-P, P-S, S-P, and S-S components, are calculated on a sparse grid using elastic VTI wave equation Born modeling and migration. The traditional migration result is then decomposed into individual local results, each corresponding to the size of a single PSF, through unit partitioning. Finally, a set of filters is derived from four elastic VTI PSF components to approximate the Hessian inverse, which are subsequently assembled to reconstruct an entire deconvolved image. The numerical experiments confirm that our method can significantly decrease multiparameter crosstalk artifacts and improve spatial resolution and amplitude fidelity.

Efficient preconditioned least-squares wave-equation migration

Since the appearance of wave-equation migration (WEM), many have tried to improve the resolution and effectiveness of this technology. Least-squares wave-equation migration is one of those attempts that try to fill the gap between migration assumptions and reality in an iterative manner. However, these iterations do not come cheap. A proven solution to limit the number of least-squares iterations is to correct the gradient direction within each iteration via the action of a preconditioner that approximates the inverse Hessian. However, the Hessian computation, or even the Hessian approximation computation, in large-scale seismic imaging problems involves an expensive computational bottleneck, making it unfeasible. Therefore, we develop an efficient computation of the Hessian approximation operator in the context of one-way WEM in the space-frequency domain. We build the Hessian approximation operator depth by depth, considerably reducing the operator size each time it is calculated. We prove the validity of our method with two numerical examples. We then extend our proposal to the framework of full-wavefield migration, which is based on WEM principles but includes interbed multiples. Finally, this efficient preconditioned least-squares full-wavefield migration is successfully applied to a data set with strong interbed multiple scattering.

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.

Full-waveform inversion by model extension: Theory, design, and optimization

We describe a new method, full-waveform inversion by model extension (FWIME), that recovers accurate acoustic subsurface velocity models from seismic data when conventional methods fail. We leverage the advantageous convergence properties of wave equation migration velocity analysis (WEMVA) with the accuracy and high-resolution nature of acoustic full-waveform inversion (FWI) by combining them into a robust mathematically consistent workflow with minimal need for user inputs. The novelty of FWIME resides in the design of a new cost function and a novel optimization strategy to combine the two techniques, making our approach more efficient and powerful than applying them sequentially. We observe that FWIME mitigates the need for accurate initial models and low-frequency long-offset data, which can be challenging to acquire. Our new objective function contains two components. First, we modify the forward mapping of the FWI problem by adding a data-correcting term computed with an extended demigration operator, whose goal is to ensure phase matching between predicted and observed data, even when the initial model is inaccurate. The second component, which is a modified WEMVA cost function, allows us to progressively remove the contributions of the data-correcting term throughout the inversion process. The coupling between the two components is handled by the variable projection method, which reduces the number of adjustable hyperparameters, thereby making our solution simple to use. We devise a model-space multiscale optimization scheme by reparameterizing the velocity model on spline grids to control the resolution of the model updates. We generate three cycle-skipped 2D synthetic data sets, each containing only one type of wave (transmitted, reflected, or refracted), and we analyze how FWIME successfully recovers accurate solutions with the same procedure for all three cases. In a second paper, we apply FWIME to challenging realistic examples where we simultaneously invert all wave modes.

Full-waveform inversion by model extension: Practical applications

Producing reliable acoustic subsurface velocity models still remains the main bottleneck of the oil and gas industry’s traditional imaging sequence. In complex geologic settings, the output of conventional ray-based or wave-equation-based tomographic methods may not be sufficiently accurate for full-waveform inversion (FWI) to converge to a geologically satisfactory earth model. We create a new method referred to as full-waveform inversion by model extension (FWIME) in which a wave-equation migration velocity analysis (WEMVA) technique is efficiently paired with a modified version of FWI. We find that our method is more powerful than applying WEMVA and FWI sequentially, and that it can converge to accurate solutions without the use of a good initial guess or low-frequency energy. We determine FWIME’s potential on five realistic and challenging numerical examples that simulate complex geologic scenarios often encountered in hydrocarbon exploration. We guide the reader step by step throughout the optimization process. We find that our method can simultaneously invert all wave types with the same simple mechanism and without the need for a user-intensive hyperparameter tuning process. In an open-source online repository, we provide our C++/compute unified device architecture (CUDA) numerical implementation accelerated with graphics processing units, encapsulated in a Python interface. All the numerical examples developed are accessible through Python notebooks and fully reproducible.

A new data-domain reflection full-waveform inversion method based on common-image-point gathers

The migration velocity analysis (MVA) based on the common-image gathers (CIGs) in the image domain is an effective approach for evaluating the seismic background velocity, a crucial index for the full-waveform inversion (FWI) and migration. Compared with cumbersome conventional MVA methods, differential semblance optimization (DSO) is a variation of image-domain wavefield tomography that avoids interactive manual picking but increases the computation cost tremendously. In recent years, the reflection waveform inversions have been studied in depth to construct the background velocity, matching reflection travel time or waveform in the data domain, which requires expensive Born simulations based on accurate imaging results. This study performs a new reflection full-waveform inversion method in the image domain based on the flatness of arbitrarily selected common-image-point (CIP) gathers and their Born simulations. According to the relationship between imaging and background model, the flatness of CIP gathers under the L1 norm is measured. Afterwards, we construct the objective function to update the background model without manual picking or wave-equation migration in an extended dimension. The inversion accuracy is controlled by balancing the density of the imaging points and the computation capability. The whole procedure is implemented in the frequency domain, so the main computation comes from the LU decomposition in the frequency-domain wave-equation simulation, in which way, the proposed data-domain reflection full-waveform inversion achieves complete automation with lower costs. The objective function analysis shows a better convex, proving the weak dependence on the starting model of this method. Besides, synthetic experiments verify the validity and the field data test demonstrates its effectiveness.

Fast hyperparameter-free spectral approach for 2D seismic data reconstruction

Reconstruction of missing seismic data is a critical procedure for subsequent applications like multiple wave suppression, wave-equation migration imaging and so on. In this paper, a fast, hyperparameter-free and sparse iterative spectral estimation approach is proposed for the reconstruction of two-dimensional seismic data of randomly missing traces. The proposed approach is based on the harmonic structure of the frequency slice of seismic data and the weighted covariance fitting criterion. Specifically, the method first iteratively estimates the spectrum of the frequency slice by solving a weighted covariance fitting problem. Then, the missing data is reconstructed by using the estimated spectrum and a linear minimum mean-squared error estimator. However, the spectral estimation depends on matrix-vector multiplications for each iteration, which has a high computational cost when the data increase to a large size. To solve this problem, a fast iterative technology is proposed by using an inverse fast Fourier transform, which fully exploits the Hermitian–Toeplitz structure of the covariance matrix and the exponential form of the steering vector and it significantly reduces the computational complexity. The proposed algorithm is hyperparameter-free, can provide super spectral resolution, and thus obtain better reconstruction performance. The experimental results of synthetic and real seismic data show that the proposed algorithm has higher reconstruction accuracy and lower computational complexity compared to other commonly used reconstruction algorithms.

Introduction to a Two‐Way Beam Wave Method and Its Applications in Seismic Imaging

AbstractRay theory gives a high‐frequency asymptotic solution of the wave equation, which has been widely used in the seismic study because of its high computational efficiency and good physical insights for elementary waves. However, the fundamental high‐frequency assumption makes it difficult to accurately simulate the finite‐frequency effect of wave propagations. On the other hand, the Gaussian beam, as one of the advanced ray‐based methods, overcomes the amplitude singularity problem near the caustics by introducing a complex‐valued paraxial approximation. But Gaussian beams only use the velocity and up to second‐order velocity derivative of central rays and thus does not accurately simulate the response of heterogeneity away from the rays. To bridge the gap between the ray theory and wave‐equation methods, we present a two‐way beam wave method and apply it to seismic imaging. A fan of central rays are first shot to extract local ray‐centered model parameters in the global Cartesian coordinates. Then, a finite‐difference method is used to solve the acoustic wave equation in the ray‐centered coordinates to simulate wave propagations near the central rays. The beam wave method can accurately simulate the response of the velocity heterogeneities in the ray tubes and therefore produces more accurate wavefields. In addition, we discretize the ray‐centered wave equation onto a non‐orthogonal grid, which considerably reduces the computational cost. We apply the two‐way beam wave method in seismic imaging and develop a beam wave reverse‐time migration. It inherits the advantages of both ray‐based and wave equation migrations, and produces clear images for complicated structures with fewer artifacts than traditional imaging methods.

Efficient wavefield separation by reformulation of two-way wave-equation depth-extrapolation scheme

In two-way wave-equation depth migration by either time extrapolation or depth extrapolation, an efficient wavefield separation scheme (to separate two-way wavefields into their one-way counterparts: down- and upgoing waves) is typically required for high-quality imaging. For a depth-extrapolation scheme, this wavefield separation is conventionally accomplished by using two wavefields, the pressure p and its derivative [Formula: see text]. However, that wavefield separation scheme presents a significant challenge in the accuracy and the computational cost for two-way wave-equation-based depth migration. To address this challenge, we have developed a new wavefield separation algorithm for a two-way wave-equation depth-extrapolation scheme, in which a new pressure wavefield q is introduced to replace the pressure derivative wavefield [Formula: see text] in boundary conditions and propagator of the conventional formulation. Because the pressure wavefield p represents the sum of down- and upgoing waves, whereas the pressure wavefield q denotes the difference between down- and upgoing waves, wavefield separation in the new algorithm becomes very simple and efficient — a summation and subtraction operation only at each depth step. The performance of our algorithm for wavefield separation is verified and demonstrated by numerical experiments conducted on several synthetic models. These experiments include examples of wavefield separation of the source and the receiver, migration tests on the 2D SEG/EAGE salt model (regarding one-way versus two-way depth extrapolation and comparison of evanescent-wave-suppression methods), and turning wave migration tests on three complex velocity models: the step model, the salt dome model, and the BP salt model.

Wave-equation migration velocity analysis via the optimal-transport-based objective function

Wave-equation migration velocity analysis builds a kinematically accurate macrovelocity model containing the large-scale structures of the model for seismic imaging or full-waveform inversion (FWI). Differential semblance optimization formulates the misfit function in the subsurface domain by applying a penalty to the unfocused subsurface-offset gathers or the nonflattened angle gathers. Such a penalty leads to gradients with strong artifacts and spurious oscillations, which then leads to slow convergence. We formulate the misfit function using the optimal transport (OT) among neighboring traces in angle-domain common-image gathers (CIGs). Specifically, we measure the unflatness of the gathers by comparing the Wasserstein distance of adjacent traces. Because CIGs are ideal to form a probability distribution, which is required by OT, it is natural to use OT to measure the time shifts between two distributions. Our objective function exhibits well-behaved convex and unimodal properties in a simple model with a single interface, and it is expected to be minimum for the correct velocity when the angle gathers are flat. A synthetic data set example on the Marmousi model indicates the ability of our method to provide an initial model for FWI that allows FWI to converge. We also apply the approach on a marine field data set to further demonstrate its effectiveness in estimating the macrovelocity model.

Two dimensional dynamically focused beam migration in weakly anisotropic media

Anisotropy exists widely in the earth medium. For anisotropic medium, the traditional isotropic migration will lead to the inaccurate images and unfocused energy. Especially for a deep target exploration, ignoring anisotropy will further affect the subsequent seismic data interpretation. As an improved ray method, Gaussian beam migration overcomes the imaging problems of caustics and shadow zones in Kirchhoff migration, and avoids the time-consuming problems of wave equation migrations. However, the imaging quality of the Gaussian beam migration is controlled by the initial beam width at the surface. In this paper, we extend Nowack (2011)’s algorithm to anisotropic media and present an anisotropic dynamically focused beam migration by modifying the propagator of the Gaussian beams. We use anisotropic kinematic and dynamic ray tracing with optimized coefficients to obtain travel times, trajectories and dynamic information. We then use dynamically focused beams for weakly anisotropic media to calculate the Green's function and the Claerbout imaging condition to obtain images. This strategy enables us to improve the imaging quality in the middle and deep layers without reducing the imaging quality in the shallow layers, which will help to overcome the limitation of the initial beam width in Gaussian beam migration. The results for a VTI fault model and the Shengli complex structure model show the accuracy and validity of the proposed method in this paper.

Wavefield separation using irreversible-migration filtering

Wavefield separation is important for eliminating unwanted components of seismic data while retaining preferred components therein; however, there have been difficulties with using traditional methods to select a proper muting window in the transformed domain. A narrow muting window can separate different components well, but it causes visible artifacts due to sudden truncations (e.g., the Gibbs phenomenon), whereas a wide muting window would, on the other hand, leave unwanted components. We have adopted separating the wavefield based on its irreversibility after migration and demigration. One-way wave-equation migration would automatically damp high-slope components in the evanescent region while accurately handling the other components. Thus, the migration velocity can be tuned to naturally trap a given range of high-slope components into the evanescent region in the migrated domain, which would be irreversible after demigration. In contrast, the other components (i.e., those outside the evanescent region) can be recovered after demigration. Our method achieves perfect muting by taking advantage of migration and demigration, and thus it avoids the manual operation of setting a muting window. As a result, our method is free of muting artifacts. We conduct numerical experiments with synthetic and field data, and the results verify the excellent performance of our method for several different kinds of wavefield separation, including linear event separation, structural noise elimination, diffraction-reflection separation, and vertical seismic profile wavefield separation. Our method integrates noise reduction and the wavefield separation, and thus it can reduce the computational cost of using different data processing schemes and avoid the related potential error accumulation.

Joint inversion of the reflectivity and the velocity model

During linearized waveform inversion (LWI), the presence of small inaccuracies in the background subsurface model can lead to unfocused seismic events in the final image. The effect on the amplitude can mislead the interpretation. We have developed a joint inversion scheme in the model domain of the reflectivity and the background velocity model. The idea is to unify the inversion of the background and the reflectivity model into a single framework instead of treating them as decoupled problems. We find that, with this method, we can obtain a better estimate of the reflectivity than that obtained with conventional LWI. Conversely, the background model is improved by the joint inversion with the reflectivity in comparison with wave-equation migration velocity analysis. We perform tests on 2D synthetics and 3D field data that demonstrate both benefits.

Extending subsurface imaging beyond ocean-bottom node coverage using multiples: A case study in deepwater Nigeria

We have evaluated the results of a receiver decimation study in a deepwater context using separated wavefield imaging (SWIM) algorithms to provide extended illumination for imaging without ocean-bottom node (OBN) positioning constraints. We carried out subsurface imaging using the SWIM imaging technique with a reduced OBN layout, and we compared the results with those from conventional one-way wave-equation migration. We found from the results that the SWIM algorithm makes it possible to reduce the OBN layout while obtaining a similar subsurface image with the same shot geometry, which allows a reduced receiver acquisition effort, offers more geometry flexibility without affecting the image quality, with a potentially significant reduction of acquisition cost and 4D processing turnaround time.

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.

Deep learning for velocity model building with common-image gather volumes

SUMMARY Subsurface velocity model building is a crucial step for seismic imaging. It is a challenging problem for conventional methods such as full-waveform inversion (FWI) and wave equation migration velocity analysis (WEMVA), due to the highly nonlinear relationship between subsurface velocity values and seismic responses. In addition, traditional FWI and WEMVA methods are often computationally expensive. In this paper, we propose to apply a deep learning technique to construct subsurface velocity models automatically from common-image gather (CIG) volumes. In our method, pairs of synthetic velocity models and CIG volumes are generated to train a convolutional neural network. Our proposed network achieves promising results on different synthetic data sets. The training performance of several commonly used loss functions is also studied.

Differential semblance optimisation based on the adaptive quadratic Wasserstein distance

Abstract As the robustness for the wave equation-based inversion methods, wave equation migration velocity analysis (WEMVA) is stable for overcoming the multipathing problem and has become popular in recent years. As a rapidly developed method, differential semblance optimisation (DSO) is convenient to implement and can automatically detect the moveout existing in common image gathers (CIGs). However, by implementing in the image domain with the target of minimising moveouts and improving coherence of the CIGs, the DSO method often suffers from imaging artefacts caused by uneven illumination and irregular observation geometry, which may produce poor velocity updates with artefact contamination. To deal with this issue, in this paper, by introducing Wiener-like filters, we modify the conventional image matching-based objective function to a new one by introducing the quadratic Wasserstein metric technique. The new misfit function measures the distance of two distributions obtained by the convolutional filters and target functions. With the new misfit function, the adjoint sources and the corresponding gradients are improved. We apply the new method to two numerical examples and one field dataset. The corresponding results indicate that the new method is robust to compensate low frequency components of velocity models.

Case study: Improving the quality of the seismic reflection image for a Fujian mineral exploration data set with offset-domain common-image gathers

Seismic reflection is a proven and effective method commonly used during exploration of deep mineral deposits in Fujian, China. In seismic data processing, rugged depth migration based on wave-equation migration can play a key role in handling surface fluctuations and complex underground structures. Because wave-equation migration in the shot domain cannot output offset-domain common-image gathers in a straightforward way, the use of traditional tools for updating the velocity model and improving image quality can be quite challenging. To overcome this problem, we have used the attribute migration method. This works by sorting the migrated stack results for every single-shot gather into the offset gathers. The value of the offset that corresponds to each image point is obtained from the ratio of the original migration results to the offset-modulated shot-data migration results. A Gaussian function is proposed to map every image point to a certain range of offsets. This helps improve the signal-to-noise ratio, which is especially important in handling low-quality seismic data obtained during mineral exploration. Residual velocity analysis is applied to these gathers to update the velocity model and improve the image quality. The offset-domain common-image gathers are also used directly for real mineral exploration seismic data with rugged depth migration. After several iterations of migration and updating the velocity, our procedure achieves an image quality better than the one obtained with the initial velocity model. The results can help with the interpretation of thrust faults and deep deposit exploration.

A multi-level parallel algorithm for seismic imaging based on one-way wave equation migration

This article presents a parallel algorithm for seismic imaging based on depth wavefield extrapolation (the solution of a one-way wave equation [OWE]) employing the pseudospectral method. The algorithm is essentially oriented on common-offset vector (COV) gathers seismic migration based on an OWE, includes several parallelism levels, and utilizes message passing interface (MPI), Nested OpenMP, and CUDA technologies. The uppermost level of parallelism involves data splitting to ensure that each dataset can be processed independently to construct parts of the COV images. The algorithm is embarrassingly parallel at this level. Next, each independent run processes several datasets to compute COV images. Each MPI process computes all COV images for a single dataset, then MPI processes exchange data to build up a single COV image for all datasets at a single node. Computation of the COV image at a node requires OpenMP parallelization so that one thread governs the GPU-based calculations and facilitates the construction of images within a thick slab. All other threads perform image interpolation within the slab. Finally, CUDA technology is used for the most computationally intense part of the algorithm—wavefield extrapolation. Such a complex algorithm structure makes it possible to process all COV images for full-azimuth seismic data employing OWE-based migration.