Stéphane Operto

SUMMARY: We present a novel model-based seismic data redatuming method based upon time-domain data-assimilated wavefield (DAW) reconstruction. Seismic redatuming refers to the positioning of sources and/or receivers from the acquisition surface to a virtual datum level along with the computation of the data set that would have been recorded by the virtual acquisition. Physically, DAWs approximate the true wavefields by least-squares time-reversal extrapolation of the recorded data in a background subsurface model. To this end, DAWs are first formulated as the sum of the background wavefields and an approximation of the scattered wavefields by the unknown contrasts between the true model and the background model. We estimate the scattered fields via the estimation of their volume scattering sources by fitting in a least-squares sense the differences between the recorded data and the simulated data in the background model, namely the restriction at receivers of the scattered fields. These underdetermined linear source problems involve two main steps: a multidimensional deconvolution of the scattered data by the data-domain Hessian, followed by the propagation of the deconvolved scattered data backward in time. Once the source problems have been solved, DAWs follow by solving the wave equation in the background model using the sum of the experimental sources and the scattering sources as extended sources. Finally, the redatumed data sets, whose accuracy depends on the precision of the background medium between the acquisition surface and the datum, are readily obtained by sampling the DAWs at arbitrary datum levels. The primary challenge in computing accurate DAWs lies in the multidimensional deconvolution of the scattered data, which requires solving a data-domain normal equation with pre-conditioned Krylov-subspace iterative methods, where each iteration requires one forward and one backward simulation. We assess the accuracy of the method and discuss its computational cost with several benchmarks representing various experimental environments (onshore with weathered layers, shallow and deep offshore).

Evolving from acoustic full-waveform inversion (AFWI) to elastic anisotropic FWI of long-offset, wide-azimuth multicomponent ocean-bottom node data is one of the most topical issues in exploration geophysics. The goal is to tend toward reservoir characterization by accounting for elastic effects in FWI. In this framework, two methodological issues are the computational overhead generated by elastic wave modeling and the management of the four components during multiparameter FWI. Although the finite-difference method on regular Cartesian grids is widespread in AFWI, this choice may be questioned in elastic FWI (EFWI), where complex bathymetry and low shear waves near the sea bottom may require unstructured meshing of the medium to reach the best trade-off between accuracy and computational cost. In this study, we perform a first assessment of the spectral element method (SEM) on unstructured hexahedral meshes as the forward engine of time-domain EFWI with a preliminary application on a low-frequency narrow-azimuth subset of the 2015 Gorgon OBN data set, offshore Australia. We show that EFWI can be applied to the pressure component to update both compressional and shear wave speeds (VP and VS) and that multiparameter inversion leads to better convergence than monoparameter inversion for VP. The good data fit achieved after FWI, the structural consistency between the long-wavelength EFWI model and the 25 Hz AFWI model reconstructed during a previous study, and the consistent VP/VS ratio in the reservoir with empirical relationships relating VP and VS for sand and shale lithologies provide a first validation of our SEM-based EFWI workflow. The next step is the inversion of the four components for joint reconstruction of the density and the full elasticity tensor.

Extended-source full-waveform inversion (ES-FWI) first computes wavefields with data-driven source extensions such that the simulated data in inaccurate velocity models match the observed counterpart sufficiently well to prevent cycle skipping. Then, the source extensions are minimized to update the model parameters toward the true medium. This two-step workflow is iterated until data and sources are matched. It has been recently indicated that the source extensions are the least-squares solutions of the scattered data-fitting problem. As a result, the source extensions are computed by propagating backward in time the deconvolved data residuals by the damped data-domain Hessian of the scattered data-fitting problem. Estimating these weighted data residuals is the main computational bottleneck of time-domain ES-FWI. To mitigate this burden, we approximate the inverse data-domain Hessian by mono- and multidimensional matching filters with two simulations per source. We implement time-domain ES-FWI with the alternating-direction method of multipliers and total-variation regularization. Moreover, we apply ES-FWI with a multiscale approach involving frequency continuation and layer stripping, with the latter being implemented with an offset-time-dependent weighting operator. In this framework, we further regularize the inversions while mitigating their computational burden by matching the grid interval to the frequency bandwidth. Finally, the overall workflow combines ES-FWI and classical FWI during the early and late stages of the multiscale approach, respectively. We illustrate that the sensitivity of ES-FWI to the accuracy of the approximated inverse data-domain Hessian depends on the complexity of the targeted model, the data anatomy, and the accuracy of the starting model. In the case of the 2004 BP salt model, we determine that the layer stripping is necessary when the inverse data-domain Hessian is approximated by a 2D Gabor matching filter and the starting model is crude, whereas this feature is not necessary with the Marmousi II model.

Full-waveform inversion (FWI) can be made immune to cycle skipping by matching the recorded data with traveltime errors smaller than one-half period from inaccurate subsurface models. To achieve this goal, the simulated wavefields can be computed in an extended search space as the solution of an overdetermined problem aimed at jointly satisfying the wave equation and fitting the data in a least-squares sense. This leads to data-assimilated wavefields that are computed by solving the wave equation in the inaccurate background model with a data-dependent source extension added to the source term. Then, the subsurface parameters are updated by canceling out these additional source terms, sometimes inaccurately called wave equation errors, to push the background model toward the true model in the left-side wave equation operator. Although many studies are devoted to these approaches with promising numerical results, their governing physical principles and their relationships with classical FWI do not seem to be understood well yet. The goal of this tutorial is to review these principles in the framework of inverse scattering theory whose governing forward equation is the Lippmann-Schwinger equation. From this equation, we find how the data-assimilated wavefields embed an approximation of the scattered field generated by the sought model perturbation and how they modify the sensitivity kernel of classical FWI beyond the Born approximation. We also clarify how the approximation with which these wavefields approximate the unknown true wavefields is accounted for in the adjoint source and in the full Newton Hessian of the parameter-estimation problem. The theory is finally illustrated with numerical examples. Understanding the physical principles governing these methods is a necessary prerequisite to assess their potential and limits and design relevant heuristics to manage the latter.

Full-waveform inversion (FWI) with extended sources first computes wavefields with data-driven source extensions, such that the simulated data in inaccurate velocity models match the observed counterpart well enough to prevent cycle skipping. Then, the source extensions are minimized to update the model parameters. This two-step workflow is iterated until both data and sources are matched. It was recently shown that the source extensions are the least-squares solutions of the recorded scattered data fitting problem. As a result, they are computed by propagating backward in time the deblurred FWI data residuals, where the deblurring operator is the inverse of the damped data-domain Hessian of the scattering-source estimation problem. Estimating the deblurred data residuals is the main computational bottleneck of time-domain extended-source FWI (ES-FWI). To mitigate this issue, we first estimate them when the inverse of the data-domain Hessians is approximated by matching filters in Fourier and short-time Fourier domains. Second, we refine them with conjugate-gradient iterations when necessary. Computing the matching filters and performing one conjugate-gradient iteration each require two simulations per source. Therefore, it is critical to design some workflows that minimize this computational burden. We implement time-domain ES-FWI with the augmented Lagrangian method. Moreover, we further extend its linear regime with a multiscale frequency continuation approach, which is combined with grid coarsening to mitigate the computational burden and regularize the inversion. Finally, we use total-variation regularization to deal with large-contrast reconstruction. We present synthetic cases where different inversion workflows carried out with data-domain Hessians of variable accuracy were assessed with the aim at converging toward accurate solutions while minimizing computational cost.

Frequency-domain full-waveform inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with a sparse layout of ocean-bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing the solution of the forward (boundary-value) problem efficiently in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with a subset of the 2015/2016 Gorgon OBN data set in the North West Shelf, Australia. We solve the forward problem with the massively parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: an efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed-precision arithmetic, and efficient processing of sparse sources. The Gorgon subdata set involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Monoparameter FWI for vertical wavespeed is performed in the viscoacoustic vertically transverse isotropic approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency ≥ 8.5 Hz) to 705 km2 (frequency ≤ 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intrahorst faults and several horizons of the Mungaroo Formation down to a depth of 7 km. Seismic modeling reveals a good kinematic agreement between recorded and simulated data, but amplitude mismatches between the recorded and simulated reflection from the reservoir suggest elastic effects. Therefore, future works involve multiparameter reconstruction for density and attenuation before considering elastic FWI from hydrophone and geophone data.

Full-waveform inversion (FWI) is a nonlinear optimization problem that addresses the estimation of subsurface model parameters by matching the predicted to the observed seismograms. We formulate FWI as a constrained optimization problem in which the regularization term is minimized subject to the nonlinear data matching constraints. Unlike standard FWI, which solves this regularized problem with a penalty method, we use an augmented Lagrangian formulation in the framework of the method of multipliers (MM). This leads to a two-step recursive algorithm, which is called multiplier waveform inversion (MWI). First, the primal step solves an unconstrained minimization problem, which is similar to the standard FWI but with a modified data term. These modified data are obtained by adding the Lagrange multipliers to the data residuals of the classical FWI. Then, the dual step updates the Lagrange multipliers with basic gradient ascent steps. In this framework, they reduce to the running sum of the data residuals of the previous iterations. The performance of the overall algorithm is improved by considering that the MM does not require the exact solution of the primal step at each iteration. In fact, convergence is attained when the primal subproblem is solved with one (without an inner loop) iteration of a gradient-based method at each iteration of the MWI. The new algorithm can be easily implemented in existing FWI codes without computational overhead by adding the Lagrange multipliers to the data residuals at each iteration. We find with numerical examples that MWI has a faster convergence speed and much improved stability compared with FWI and can converge to an accurate solution of the inverse problem in the absence of low-frequency data, even with a constant step size.

Full Waveform Inversion can be made immune to cycle skipping by matching the recorded data arbitrarily well from inaccurate subsurface models. To achieve this goal, the simulated wavefields can be computed in an extended search space as the solution of an overdetermined problem aiming at jointly satisfying the wave equation and fitting the data in a least-squares sense. Simply put, the wavefields are computed by solving the wave equation in the inaccurate background model with a feedback term to the data added to the physical source in the right-hand side. Then, the subsurface parameters are updated by canceling out these additional source terms, sometimes called unwisely wave-equation errors, to push the background model toward the true model in the left-hand side wave-equation operator. Although many studies were devoted to these approaches with promising numerical results, their governing physical principles and their relationships with classical FWI don't seem to be understood well yet. The goal of this tutorial is to review these principles in the theoretical framework of inverse scattering theory whose governing forward equation is the Lippmann-Schwinger equation. From this equation, we show how the data-assimilated wavefields embed an approximation of the scattered field generated by the sought model perturbation and how they modify the sensitivity kernel of classical FWI beyond the Born approximation. We also clarify how the approximation with which these wavefields approximate the unknown true wavefields is accounted for in the adjoint source of the parameter estimation problem. The theory is finally illustrated with numerical examples. Understanding the physical principles governing these methods is a necessary prerequisite to assessing their potential and limits and designing relevant heuristics to manage the latter.

Frequency-domain Full Waveform Inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with sparse layout of ocean bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing efficiently the solution of the forward (boundary-value) problem in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with the 2015/16 Gorgon OBN case study in the NorthWestern shelf, Australia. We solve the forward problem with the massively-parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: An efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed precision arithmetic and efficient processing of sparse sources. The Gorgon subdataset involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Mono-parameter FWI for vertical wavespeed is performed in the visco-acoustic VTI approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency > 8.5 Hz) to 705 km2 (frequency < 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intra-horst faults and several horizons of the Mungaroo formation down to a depth of 7 km.

Frequency-domain simulation of seismic waves plays an important role in seismic inversion, but it remains challenging in large models. The recently proposed physics-informed neural network (PINN), as an effective deep learning method, has achieved successful applications in solving a wide range of partial differential equations (PDEs), and there is still room for improvement on this front. For example, PINN can lead to inaccurate solutions when PDE coefficients are non-smooth and describe structurally-complex media. In this paper, we solve the acoustic and visco-acoustic scattered-field wave equation in the frequency domain with PINN instead of the wave equation to remove source singularity. We first illustrate that non-smooth velocity models lead to inaccurate wavefields when no boundary conditions are implemented in the loss function. Then, we add the perfectly matched layer (PML) conditions in the loss function of PINN and design a quadratic neural network to overcome the detrimental effects of non-smooth models in PINN. We show that PML and quadratic neurons improve the results as well as attenuation and discuss the reason for this improvement. We also illustrate that a network trained during a wavefield simulation can be used to pre-train the neural network of another wavefield simulation after PDE-coefficient alteration and improve the convergence speed accordingly. This pre-training strategy should find application in iterative full waveform inversion (FWI) and time-lag target-oriented imaging when the model perturbation between two consecutive iterations or two consecutive experiments can be small.