An open-source framework for the implementation of large-scale integral operators with flexible, modern high-performance computing solutions: Enabling 3D Marchenko imaging by least-squares inversion

Abstract

Numerical integral operators of convolution type form the basis of most wave-equation-based methods for processing and imaging of seismic data. As several of these methods require the solution of an inverse problem, multiple forward and adjoint passes of the modeling operator are generally required to converge to a satisfactory solution. This work highlights the memory requirements and computational challenges that arise when implementing such operators on 3D seismic datasets and their usage for solving large systems of integral equations. A Python framework is presented that leverages libraries for distributed storage and computing, and provides a high-level symbolic representation of linear operators. A driving goal for our work is not only to offer a widely deployable, ready-to-use high-performance computing (HPC) framework, but to demonstrate that it enables addressing research questions that are otherwise difficult to tackle. To this end, the first example of 3D full-wavefield target-oriented imaging, which comprises of two subsequent steps of seismic redatuming, is presented. The redatumed fields are estimated by means of gradient-based inversion using the full dataset as well as spatially decimated versions of the dataset as a way to investi-gate the robustness of both inverse problems to spatial aliasing in the input dataset. Our numerical example shows that when one spatial direction is finely sampled, satisfactory redatuming and imaging can be accomplished also when the sampling in other direction is coarser than a quarter of the dominant wavelength. While aliasing introduces noise in the redatumed fields, they are less sensitive to well-known spurious artefacts compared to cheaper, adjoint-based redatuming techniques. These observations are shown to hold for a relatively simple geologic structure, and while further testing is needed for more complex scenarios, we expect them to be generally valid while possibly breaking down for extreme cases

Accelerating Seismic Redatuming Using Tile Low-Rank Approximations on NEC SX-Aurora TSUBASA

“…Note that due to the extremely large size of these matrices, whilst the problem is written in a compact matrix-vector formulation, its numerical solution is performed using matrix-free operators and iterative solvers such as conjugate gradient least-squares (CGLS) or LSQR [29]. Focusing now our attention on the multi-dimensional convolution (MDC) integral operator, which represents the most expensive computations in the overall chain of operations, its inner working can be written more explicitly as follows:…”

Time-Domain Multidimensional Deconvolution: A Physically Reliable and Stable Preconditioned Implementation

“…Though our examples in this paper are 2D, our approach is general and immediately applicable to 3D seismic data. In a companion paper, Ravasi and Vasconcelos [30] show an HPC-ready setup of the very same PyLops implementation we use here and demonstrate it with 3D MDD in the image domain after redatuming. Moreover, the generality of our approach also extends directly to the cases of MDD in elastic [11,37,39] and attenuative media [12], with the only difference being that, of course, our proposed preconditioning operators must be appropriately re-parameterized for those cases (e.g., see discussion in the next paragraph).…”

“…Such a solver only requires the implementation of the forward and adjoint operations of the modelling and preconditioning operators on a given data vector. From an implementation perspective, the multi-dimensional convolution is more efficiently solved in the Fourier domain, thus, an equivalent time domain representation requires the definition of an operator Q = F −1 QF acting on a given vector and performs a step of batch matrix multiplication as described by Ravasi and Vasconcelos [30]. Here F and F −1 , are forward and inverse space-time Fourier transform operators.…”

“…Alternatively, we resolve to set up the inversion on the time-space domain and opt for a gradient-based iteration solver (e.g., LSQR). While an explicit definition in computer memory of the discrete wavefield operator P + mko is prohibitive, such solvers only require the resulting forward and adjoint operations on the data vector; therefore, a Pylops linear operator is defined instead [30,52]. Even though no preconditioned is introduced, P = I, the synthesized virtual source-receiver response, second panel in Figure 6a,b, exhibits significantly fewer spurious events, a radiating pattern with improved causality properties, plus a much closer prediction of both phase and amplitudes.…”

“…Alternatively, MDD can also be solved in the time domain as a single coupled inverse problem: whilst this approach is less known and not widely adopted yet, it has recently been shown to be more robust as solving for the entire frequency spectrum at once acts as a natural regularizer and does not require the identification of damping parameters for each frequency independently. In this direction, pioneering work has been carried out by van der Neut and Herrmann [27], van der Neut et al [28], Luiken and van Leeuwen [29], Ravasi and Vasconcelos [30], who also proposed the use of different regularization and preconditioning strategies to further alleviate the ill-posedness of the inverse problem. More specifically, van der Neut and Herrmann [27] investigate the use of sparsifying transforms (i.e., Curvelets) in combination with sparsity-promoting inversion to improve the robustness of the solution to noise in the data and mitigates artefacts arising from limited-aperture source arrays.…”