  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 modelling operator must be performed to
converge to a satisfactory solution. This work highlights the challenges that
arise when implementing such operators on 3D seismic datasets and it provides
insights into their usage for solving large systems of integral equations. A
Python framework is presented that leverages libraries for distributed storage
and computing, and provides an high-level symbolic representation of linear
operators. To validate its effectiveness, the forward and adjoint
implementations of a multi-dimensional convolution operator are evaluated with
respect to increasing size of the kernel and number of computational resources.
Our computational framework is further shown to be suitable for both classic
on-premise High-Performance Computing and cloud computing architectures. An
example of target-oriented imaging of a 3D synthetic dataset which comprises of
two subsequent steps of seismic redatuming is finally presented. In both cases,
the redatumed fields are estimated by means of least-squares inversion using
the full dataset as well as spatially decimated versions of the dataset as a
way to investigate the robustness of both inverse problems to spatial aliasing
in the input dataset. We observe that less strict sampling requirements apply
in three dimensions for these algorithms compared to their two dimensions
counterparts. Whilst aliasing introduces noise in the redatumed fields, they
are however deprived of the well-known spurious artefacts arising from
incorrect handling of the overburden propagation in cheaper, adjoint-based
redatuming techniques.
