Least-squares Migration Research Articles

Robust full waveform inversion with deep Hessian deblurring

SUMMARY Full waveform inversion (FWI) is a technique widely used in geophysics to obtain high-resolution subsurface velocity models from waveform seismic data. Due to its large computation cost, most flavours of FWI rely only on the computation of the gradient of the loss function to estimate the update direction, therefore ignoring the contribution of the Hessian. Depending on the level of computational resources one can afford, an approximate of the inverse of the Hessian can be calculated and used to speed up the convergence of FWI towards the global (or a plausible local) minimum. In this work, we propose to use an approximate Hessian computed from a linearization of the wave equation as commonly done in least-squares migration. More precisely, we rely on the link between a migrated image and a doubly migrated image (i.e. an image obtained by demigration–migration of the migrated image) to estimate the inverse of the Hessian. However, instead of using non-stationary compact filters to link the two images and approximate the Hessian, we propose to use a deep neural network to directly learn the mapping between the FWI gradient (output) and its Hessian (blurred) counterpart (input). By doing so, the network learns to act as an approximate inverse Hessian: as such, when the trained network is applied to the FWI gradient, an enhanced update direction is obtained, which is shown to be beneficial for the convergence of FWI. The weights of the trained (deblurring) network are then transferred to the next FWI iteration to expedite convergence. We demonstrate the effectiveness of the proposed approach on two synthetic data sets and a field data set.

Seismic Imaging of the Arctic Subsea Permafrost Using a Least-Squares Reverse Time Migration Method

High-resolution seismic imaging allows for the better interpretation of subsurface geological structures. In this study, we employ least-squares reverse time migration (LSRTM) as a seismic imaging method to delineate the subsurface geological structures from the field dataset for understanding the status of Arctic subsea permafrost structures, which is pertinent to global warming issues. The subsea permafrost structures in the Arctic continental shelf, located just below the seafloor at a shallow water depth, have an abnormally high P-wave velocity. These structural conditions create internal multiples and noise in seismic data, making it challenging to perform seismic imaging and construct a seismic P-wave velocity model using conventional methods. LSRTM offers a promising approach by addressing these challenges through linearized inverse problems, aiming to achieve high-resolution, subsurface imaging by optimizing the misfit between the predicted and the observed seismic data. Synthetic experiments, encompassing various subsea permafrost structures and seismic survey configurations, were conducted to investigate the feasibility of LSRTM for imaging the Arctic subsea permafrost from the acquired seismic field dataset, and the possibility of the seismic imaging of the subsea permafrost was confirmed through these synthetic numerical experiments. Furthermore, we applied the LSRTM method to the seismic data acquired in the Canadian Beaufort Sea (CBS) and generated a seismic image depicting the subsea permafrost structures in the Arctic region.

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.

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.

High-precision seismic imaging for complex deep structures in the hydrocarbon exploration using a coherent-stacking-based least-squares migration

High-precision seismic imaging for complex deep structures in the hydrocarbon exploration using a coherent-stacking-based least-squares migration

Deep-unrolling architecture for image-domain least-squares migration

Deep-image prior (DIP) is a novel approach to solving ill-posed inverse problems whose solution is parameterized with an untrained deep neural network and cascaded with the forward modeling operator. A key component to the success of such a method is represented by the choice of the network architecture, which must act as a natural prior to the inverse problem at hand and provide a strong inductive bias toward the desired solution. Inspired by the close link between neural networks and iterative algorithms in classical optimization, we apply an unrolled version of the gradient descent (GD) algorithm as our DIP network architecture, denoted as the deep-unrolling (DU) architecture. Each layer of the unrolled network is comprised of two parts: the first part corresponds to the GD step of the data-fidelity term, whereas the second part, formed by a six-layer convolutional neural network (CNN), plays the role of a regularizer in the associated objective function. The developed DU architecture is applied to the problem of image-domain least-squares migration (IDLSM) to invert migrated seismic images for their underlying reflectivity and is denoted as DU-IDLSM. As such, the DU architecture parameterizes the reflectivity, and the input of each layer of the unrolled network is the reflectivity at the previous layer. Similar to the classical DIP approach, the parameters of the DU architecture are optimized in an unsupervised fashion by minimizing the data misfit function itself. Through experiments with a part of the Sigsbee2A model and a marine field data set, we test the effectiveness of the DU-IDLSM approach and highlight two key benefits. First, the DU architecture can effectively regularize the inversion process, resulting in reflectivity estimates with fewer artifacts and higher image resolution than those produced by conventional IDLSM approaches. Second, we indicate that by including dropout layers in the CNN architecture, DU-IDLSM can produce a qualitative measure of the uncertainty associated with the least-squares migration process.

Efficient least-squares reverse time migration in TTI media using a finite-difference solvable pure qP-wave equation

Abstract The anisotropic characteristics of underground media play a crucial role in seismic wave propagation and should be accounted for during seismic imaging. The least-squares reverse time migration (LSRTM) method achieves ideal imaging results by suppressing migration artifacts and balancing amplitude. Therefore, pure quasi-P (P-qP) wave equations in tilted transversely isotropic (TTI) media have been widely used to implement LSRTM, owing to their ability to produce stable and noise-free wavefields. However, solving the anisotropic P-qP-wave equations typically necessitates the use of spectral-based methods, making them computationally inefficient, especially in 3D applications. In this study, we first develop a P-qP-wave equation in TTI media that can be efficiently computed through the finite-difference (FD) approach. Numerical tests show that, in comparison to the previous TTI P-qP-wave equation, the newly derived FD solvable TTI P-qP-wave equation yields reasonably accurate and highly efficient wavefield simulations. Then, building on our newly derived wave equation, we derive its adjoint migration and demigration operator to implement TTI LSRTM. Two synthetic examples suggest that the newly presented P-qP-wave TTI LSRTM approach effectively correct for the anisotropy effects, providing high-quality imaging results. Additionally, our approach has superior computational efficiency over the conventional P-qP-wave TTI LSRTM technique based on a hybrid FD pseudo-spectral (HFDPS) solver.

Full-Wavefield Migration Using an Imaging Condition of Global Normalization Multi-Order Wavefields: Application to a Synthetic Dataset

In marine seismic exploration, seismic signals comprise primaries that undergo first-order scattering, as well as multiples resulting from multi-order scattering events. Surface-related multiples involve multi-order scattering at the free surface interface between seawater and air and exhibit a smaller reflection angle and broader illumination compared to primaries. Internal multiples, originating from multi-order scattering among stratified layers, provide additional illumination compensation beneath the reflecting interface. However, in conventional primary migration, different-order wavefields may result in crosstalk artifacts. To address this issue, we developed a least-squares migration (LSM) method based on the multi-order wavefield global normalization condition. This methodology investigates the illumination effects and crosstalk artifacts associated with different-order surface-related and internal multiples, and then modifies the global normalization condition by incorporating an illumination compensation perspective. Virtual sources, represented by surface-related multiples and internal multiples, are integrated into the source compensation term, ultimately yielding a multi-order wavefield normalization condition. This normalization condition is subsequently combined with least-squares full-wavefield migration (LSFWM). Numerical experiments demonstrate that the normalization condition of multi-order wavefields can resolve the problem of weak deep imaging energy and promote the suppression of multiple crosstalk artifacts in the least-squares algorithm.

Time-Lapse One-Step Least-Squares Migration

Time-Lapse One-Step Least-Squares Migration

Nonlinear least-squares reverse time migration of prismatic waves for delineating steeply dipping structures

Prismatic reflections in seismic data carry abundant information about subsurface steeply dipping structures, such as salt flanks or near-vertical faults, playing an important role in delineating these structures when effectively used. Conventional linear least-squares reverse time migration (L-LSRTM) fails to use prismatic waves due to the first-order Born approximation, resulting in a blurry image of steep interfaces. We develop a nonlinear LSRTM (NL-LSRTM) method to take advantage of prismatic waves for the detailed characterization of subsurface steeply dipping structures. Compared with current least-squares migration methods of prismatic waves, our NL-LSRTM is nonlinear and thus avoids the challenging extraction of prismatic waves or the prior knowledge of L-LSRTM results. The gradient of NL-LSRTM consists of the primary and prismatic imaging terms, which can accurately project observed primary and prismatic waves into the image domain for the simultaneous depiction of near-horizontal and near-vertical structures. However, we find that the full Hessian-based Newton normal equation has two similar terms, which prompts us to make further comparison between the Newton normal equation and our NL-LSRTM. We determine that the Newton normal equation is problematic when applied to the migration problem because the primary reflections in the seismic records will be incorrectly projected into the image along the prismatic wavepath, resulting in an artifact-contaminated image. In contrast, the nonlinear data-fitting process included in the NL-LSRTM contributes to balancing the amplitudes of primary and prismatic imaging results, thus making NL-LSRTM produce superior images compared with the Newton normal equation. Several numerical tests validate the applicability and robustness of NL-LSRTM for the delineation of steeply dipping structures and illustrate that the imaging results are much better than the conventional L-LSRTM.

Staining algorithm for least-squares reverse time migration

The staining algorithm-based reverse time migration (RTM) technique has shown promise in improving imaging results for certain complex structures, such as subsalt and steeply dipping layers, by incorporating the complex domain. However, RTM suffers from low-resolution images as it utilizes the adjoint of the forward modeling operator rather than its inverse operator. To address this limitation, we propose a new approach called stained least-squares RTM (ST-LSRTM) that enhances the gradient through the integration of the staining algorithm. In ST-LSRTM, we first develop the wave equation and Born approximation, transitioning from the real to the complex domains. To simplify the linear perturbation operator, we constrain the magnitude of the imaginary velocity. Next, we compute the gradient in the complex domain. The enhanced gradient is obtained by replacing the traditional gradient with the normalized stained gradient within the stained region. Finally, we construct a misfit function based on the L2 norm and calculate the step size and gradient direction using the conjugate gradient method. Model tests demonstrate that ST-LSRTM achieves high-resolution imaging results even when the stained area is imprecise. Comparisons of waveform curves and normalized residual curves confirm that ST-LSRTM approaches true reflection coefficient model more closely compared to traditional LSRTM.

Newton optimization and the Hessian matrix in scale-separation waveform inversion problems with reflected and transmitted waves

Summary Full waveform inversion (FWI) is often utilized for determining the medium parameters of the Earth by minimizing the misfit between observed and modelled waves. The conventional formulation of FWI is well suited to utilize transmitted waves to determine long-to-intermediate wavelength properties of the medium. Hovewer, the penetration depth of transmitted waves may be limited for a given acquisition. As a remedy, reflection waveform inversion (RWI) exploits the kinematic information in reflected waves to retrieve these long-wavelength properties, particularly below the maximum illumination depth of transmitted waves. Previous works have demonstrated the benefits of utilizing both transmitted and reflected waves, thereby considering a reconciliation of conventional FWI and RWI. One such formulation is known as joint full waveform inversion (JFWI), where reflected and transmitted waves are separately compared in the misfit functional. A disadvantage of gradient-based optimization schemes in RWI and JFWI is the presence of strong crosstalk due to mapping of reflected data residuals along the entire reflection wavepaths. Because RWI and JFWI require frequent least-squares migration, we are incentivized to determine well-focused search directions that allow large movements between the migration steps. Thus, we propose an improvement to RWI and JFWI by utilizing Newton methods to directly account for crosstalk between parameters. The crosstalk effects are described by the second-order sensitivities of the misfit functional, i.e. the Hessian matrix. We adopt the adjoint-state method to compute the first- and second order sensitivities of the misfit functional for any of these waveform inversion procedures in a unified framework and for arbitrary misfit functionals. Furthermore, we conceptually study the Hessian matrix in least-squares based full-, reflection- and joint full waveform inversion from a data-domain point of view. Synthetic examples and a three-dimensional field data example demonstrate that the truncated Gauss-Newton method is effective in attenuating artefacts due to crosstalk between spatially separated variables, thus substantially improving the spatial focusing of the search directions.

Source-independent least-squares reverse time migration in vertical transversely isotropic media based on the Student’s t-distribution

Abstract Seismic anisotropy exists in various types of strata and should be considered in seismic imaging schemes. Seismic imaging algorithms based on isotropic assumption neglect the effects of anisotropy on seismic data, which cause migration artifacts and waveform distortion. To correct the effects of anisotropy on seismic wave propagation, we propose an imaging algorithm that performs least-squares reverse time migration in vertical transversely isotropic acoustic media. We derive the following operators to implement this algorithm, the de-migration operator, its adjoint migration operator, and the corresponding gradient. However, an inaccurate estimated source wavelet will introduce error in the seismic simulation, and thus increase the mismatch between observed and synthetic data for least-squares reverse time migration. In addition, the noises, especially the noises with abnormal amplitudes in the seismic data, damage the inversion convergence and reduce the imaging resolution. To improve the image quality, we propose to use convolved wavefields between observed and synthetic data so that such a mismatch can be independent of the source wavelets. Also, we use the Student’s t-distribution instead of L2 norm in our inversion scheme to better handle the seismic noise. Its implementation only modifies the gradient of the conventional least-square reverse time migration scheme. Our numerical tests show a clear improvement using our proposed imaging algorithm when compared with the conventional isotropic migration scheme for the anisotropic data. Also, the synthetic examples demonstrate the feasibility and effectiveness of our proposed source-independent algorithm using the Student’s t-distribution.

Expectation-maximization framework-based image-domain least-squares migration

High-resolution and high-fidelity seismic inversion imaging is one of the core issues of lithologic reservoir exploration. As a linearization subproblem of full-waveform inversion, least-squares migration (LSM) can obtain inversion results with more reliable amplitude, better resolution, and fewer artifacts that can be solved in either the data domain or image domain. For data-domain LSM, the calculation of the data-domain iteration can be very computationally intensive. In practical situations, the convexity of the cost function of LSM could be poor, resulting in slow or even nonconvergence iteration, so the ideal inversion imaging results cannot be reached. The image-domain LSM (ID-LSM) follows the logic of image deblurring, which realizes the effects of LSM in another way. The key problem of ID-LSM is to calculate a reliable Hessian matrix that can be effectively replaced by a set of point-spread functions (PSFs). Affected by the unknown source wavelet, the accuracy of the background velocity, and the accuracy of the forward/migration operator, the calculation results of the Hessian/PSF cannot be accurate. We first illustrate the impact of incorrect Hessian/PSF on ID-LSM through numerical experiments. To eliminate the influence of the inaccuracy of PSF, we develop the ID-LSM problem as a high-dimensional image deconvolution problem. An alternate update strategy of PSF and imaging results is given under the framework of the expectation-maximization algorithm by additionally introducing prior information about PSF and the inversion result. Numerical examples of synthetic data and field data verify that our method significantly improves image quality.

Plane-wave least-squares diffraction imaging using short-time singular spectrum analysis

Abstract Diffractions are seismic waves generated by small-scale heterogeneities in the subsurface. These are often superimposed by strong reflections so that they are not visible on the image, leading to misinterpretation and incorrect localization of the scatterers. Therefore, the separation of diffracted and reflected waves is a crucial step in identifying these small-scale diffractors. To realize the separation of diffraction and imaging, a least-squares reverse time migration method of plane waves (PLSRTM) optimized with short-time singular spectrum analysis (STSSA) was developed in this work. The proposed STSSA algorithm exploits the properties of singular spectral analysis (SSA) to separate linear signals. By establishing the Hanning window and the energy compensation function, it also compensates for the shortcomings of SSA in local dip processing and convergence of linear signals. As there is no clear boundary between reflected and diffracted waves, the energy loss during separation leads to a slow convergence rate of the diffraction wave imaging technique. We use STSSA as a constraint for PLSRTM, which greatly improves the imaging quality for diffraction waves. The tests with the Sigsbee2A model and noisy seismic data have shown that our method can effectively improve the resolution of diffraction wave imaging and that the constraint of STSSA increases the robustness to noisy data.

Target-oriented elastic full-waveform inversion through acoustic extended image-space redatuming

Elastic full-waveform inversion (FWI) can provide accurate and high-resolution subsurface parameters. However, its high computational cost prevents the application of this method to large-scale field-data scenarios. To mitigate this limitation, we have developed a target-oriented elastic FWI methodology based on a redatuming step that relies upon an extended least-squares migration process. In our approach, the surface-reflection data can be attributed to a given subsurface portion when mapped into the image space. This process allows us to reconstruct reflection data generated by a target area and recorded with a virtual acquisition geometry positioned directly above it. The redatuming step enables the application of an elastic FWI method within the target portion only. The entire workflow drastically decreases the overall cost of the surface-data inversion and allows the retrieval of accurate elastic parameters of the area of interest. We determine the effectiveness of our approach on a synthetic case based on the well-known Marmousi2 model and on 3D ocean-bottom node pressure data recorded in the Gulf of Mexico in which we retrieve the elastic parameters of a potential prospect, positioned in proximity of a salt-dome flank, and whose rock-physical properties are consistent with the presence of a gas-bearing sand reservoir.

Improved least-squares migration through double-sweeping solver

Based on a recently developed approximate wave-equation solver, we have developed a methodology to reduce the computational cost of seismic migration in the frequency domain. This approach divides the domain of interest into smaller subdomains, and the wavefield is computed using a sequential process to determine the downward- and upward-propagating wavefields — hence called a double-sweeping solver. A sequential process becomes possible using a special approximation of the interface conditions between subdomains. This method is incorporated into the least-squares migration framework as an approximate solver. The associated computational effort is comparable to one-way wave-equation approaches, yet, as illustrated by the numerical examples, the accuracy and convergence behavior are comparable to that of the full-wave equation.

Phase-encoding-based least-squares reverse time migration of controlled-order multiples

Reverse time migration of controlled-order multiples (RTM-CM) can significantly remove strong crosstalk artifacts in multiples imaging. However, the computational cost is demanding if all pairs of consecutive-order multiples are considered, and an RTM-CM image retains residual crosstalk. To address these issues, we have developed a phase-encoding-based migration approach. This approach first phase-encodes different-order multiples by randomly modifying time shifts and polarity reversals, and then blends the encoded data into supergathers. As a result of using supergathers, migrations of different-order multiples are simultaneously accomplished, which considerably improves computational efficiency. As the simultaneous migration inherits the advantages of RTM-CM and provides good gradients for least-squares migration (LSM), we introduced an iterative LSM to gradually attenuate residual crosstalks in our approach. Numerical experiments on synthetic and field data sets revealed that our approach can significantly enhance imaging quality by suppressing crosstalks, increasing spatial resolutions, and reducing computational costs.

3D least-squares reverse time migration in VTI media based on pseudoacoustic wave equation and multi-GPU parallel acceleration

As the high-density seismic acquisition is progressively becoming an essential mode of seismic exploration, developing high-precision imaging methods for 3D datasets has become an urgent industry requirement. Unfortunately, the two-dimensional least-squares reverse time migration (LSRTM) can hardly deal with the influence of the third dimension because the structure of reflectors is constantly changing in nearly all directions. Furthermore, anisotropy is very common in underground media, which makes the LSRTM based on isotropic assumption fail to locate the reflectors accurately and seriously limits the identification of geological stratification. Considering the above two problems, this paper derives the least-squares reverse time migration in three-dimensional vertical transversely isotropic media (3D-VTI-LSRTM) and develops a workflow which can realize 3D-LSRTM in VTI media based on this theory. Our workflow uses multi-GPU acceleration and heterogeneous CPU/GPU parallelism. Compared with the general parallel CPU workflow or single GPU workflow, it is feasible that the GPU application significantly improves computing efficiency in practical applications. In addition, this paper extends VTI-LSRTM from 2D to 3D, which solves the problem that the 2D algorithm cannot eliminate the lateral direction reflection artifacts. For the self-made 3D-VTI depression model and the 3D-VTI-SEG/EAGE salt model, we use four NVIDIA Tesla k40c graphics cards for testing. The results of the tests are correct and significantly improved compared to the 3D-VTI-RTM method.

Building adjoint operators for least-squares migration using the acoustic wave equation

Seismic imaging can be solved as an inversion problem, which can be implemented as a least-squares migration (LSM). Compared with conventional migration algorithms, an LSM can produce imaging results with enhanced illumination and resolution. However, solving an inversion problem faces difficulties in convergence, stability, and computational efficiency. To address these issues, efforts have been spent on examining the key elements in an LSM, including the modeling operators, the migration operators, the inversion solvers, etc. Advanced modeling operators are developed to accurately simulate the seismic data. Innovative migration algorithms are implemented to precisely compute the subsurface image. Fast inversion solvers are used to improve computation efficiency. Although all these techniques are robust to improve the LSMs, they often are studied independently. It is not often considered whether these elements are properly combined when used in an LSM. We have constructed LSMs with adjoint modeling and migration operators, and we develop algorithms to prepare input shot data for these LSMs.