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 the convolution type form the basis of most wave-equation-based methods for the processing and imaging of seismic data. Because 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. We have highlighted the memory requirements and computational challenges that arise when implementing such operators on 3D seismic data sets and their use for solving large systems of integral equations. A Python framework is presented that leverages libraries for distributed storage and computing, and it 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 framework, but also 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 two subsequent steps of seismic redatuming, is presented. The redatumed fields are estimated by means of gradient-based inversion using the full data set as well as spatially decimated versions of the data set as a way to investigate the robustness of both inverse problems to spatial aliasing in the input data set. Our numerical example shows that when one spatial direction is finely sampled, satisfactory redatuming and imaging can also be accomplished when the sampling in the other direction is coarser than a quarter of the dominant wavelength. Although aliasing introduces noise in the redatumed fields, they are less sensitive to well-known spurious artifacts compared to cheaper, adjoint-based redatuming techniques. These observations are shown to hold for a relatively simple geologic structure, and although further testing is needed for more complex scenarios, we expect them to be generally valid while possibly breaking down for extreme cases (e.g., highly heterogeneous/scattering media).