PINNslope: seismic data interpolation and local slope estimation with
  physics informed neural networks
Interpolation of aliased seismic data constitutes a key step in a seismic
processing workflow to obtain high quality velocity models and seismic images.
Building on the idea of describing seismic wavefields as a superposition of
local plane waves, we propose to interpolate seismic data by utilizing a
physics informed neural network (PINN). In the proposed framework, two
feed-forward neural networks are jointly trained using the local plane wave
differential equation as well as the available data as two terms in the
objective function: a primary network assisted by positional encoding is tasked
with reconstructing the seismic data, whilst an auxiliary, smaller network
estimates the associated local slopes. Results on synthetic and field data
validate the effectiveness of the proposed method in handling aliased (coarsely
sampled) data and data with large gaps. Our method compares favorably against a
classic least-squares inversion approach regularized by the local plane-wave
equation as well as a PINN-based approach with a single network and
pre-computed local slopes. We find that introducing a second network to
estimate the local slopes whilst at the same time interpolating the aliased
data enhances the overall reconstruction capabilities and convergence behavior
of the primary network. Moreover, an additional positional encoding layer
embedded as the first layer of the wavefield network confers to the network the
ability to converge faster improving the accuracy of the data term.
Francesco Brandolin, Matteo Ravasi, Tariq Alkhalifah
PINNslope: seismic data interpolation and local slope estimation with physics informed neural networks

###### Abstract

Interpolation of aliased seismic data constitutes a key step in a seismic processing workflow to obtain high quality velocity models and seismic images. Leveraging on the idea of describing seismic wavefields as a superposition of local plane waves, we propose to interpolate seismic data by utilizing a physics informed neural network (PINN). In the proposed framework, two feed-forward neural networks are jointly trained using the local plane wave differential equation as well as the available data as two terms in the objective function: a primary network assisted by positional encoding is tasked with reconstructing the seismic data, whilst an auxiliary, smaller network estimates the associated local slopes. Results on synthetic and field data validate the effectiveness of the proposed method in handling aliased (sampled coarsely) data and data with large gaps. Our method compares favorably against a classic least-squares inversion approach regularized by the local plane-wave equation as well as a PINN-based approach with a single network and pre-computed local slopes. We find that by introducing a second network to estimate the local slopes whilst at the same time interpolating the aliased data, the overall reconstruction capabilities and convergence behavior of the primary network is enhanced. An additional positional encoding, embedded as a network layer, confers to the network the ability to converge faster improving the accuracy of the data term.

## 1 Introduction

The idea of describing seismic data as a superposition of local plane-waves was introduced by Jon Claerbout back in 1992. This elementary type of waves can be modelled by the plane-wave partial differential equation (plane-wave PDE), which is only parameterized by the local slope factor (also referred to as slowness, or ray parameter). [1] demonstrated how a plane-wave defined by its local slope can be annihilated within a given wavefield by the mean of a plane-wave partial differential operator. Leveraging on this simple concept, be developed a linear approach for local-slope estimation using small moving windows across data, where a single slope at the center of the window is computed by linear least-squares. Later, [2] proposed the plane-wave destruction filters (PWD), a global approach for local slope estimation that requires the solution of a non-linear system to estimate the dip, but removes the need for windowing the data. Several techniques have been developed during the years that utilizes the plane-wave approximation in seismic processing algorithms with applications ranging from denoising ([3]), trace interpolation, detection of local discontinuities ([2]), velocity-independent imaging ([4]) and as a regularization in seismic estimation problems ([5]).

In this work we build upon the concept introduced by Claerbout, utilizing neural networks informed by the local plane-wave equation to simultaneously interpolate seismic data and estimate the local slopes of the events. In the machine learning literature the idea of integrating the governing laws of physics into the learning process of a neural network is commonly referred to as: physics informed neural networks (PINNs - [6]). PINNs have recently emerged as a novel computational paradigm in the field of scientific machine learning and have been shown to be very effective in representing solutions of partial differential equations (PDEs). The PINNs framework can be utilized to solve both forward and inverse problems and has been successfully applied in various domains of computational physics. In the context of exploration seismology, PINNs have been utilized with different PDEs, for modelling wave propagation in time domain using the wave equation ([7], [8]) and to model wavefields in frequency domain leveraging the Helmholtz
equation and the vertical transversely isotropic wave equation ([9], [10] and [11]). PINNs have also been applied to the eikonal equation to help overcome some limitations of conventional techniques in seismic tomography problems ([12], [13], [14], [15]).

PINNs are a special class of feed-forward neural networks that predicts a solution that satisfies the PDE, which governs the physics of the problem, restricting the space of possible solutions in favor of physically reliable ones. In the standard scenario, where the loss function includes a data fitting term (as in our implementation), the PDE term acts as a soft constrain of the network optimization problem. Specifically, PINNs are generic function approximators promoting reconstructions that are naturally similar to the available data. With this advantage we can obtain an algorithm that leverages on the interpolation capabilities of feed-forward neural networks, but whose solution is heavily dependent on the physics of the problem. Moreover, compared to traditional numerical approaches that rely on the discretization of derivatives involved in a PDE, PINNs learn a direct mapping from spatial coordinates to wavefield amplitudes (or any other physical quantity), removing the need for finite difference approximations, and relying on functional derivatives, which are more stable.

A recent deep learning framework named coordinate based learning has emerged in the literature ([16], [17]). Coordinate based learning aims to solve imaging inverse problems by utilizing a feed-forward network to learn a continuous mapping from the measurement coordinates to the corresponding sensor responses. The approach bear resemblance to PINNs except that it does not utilize any physics laws to constrain the network solution. [16] found that having the feed-forward network acting directly on the input coordinates performs poorly at representing high-frequency component of the signals in accord to the low-frequency bias of neural networks as demonstrated by [18]. To overcome this problem he proposed to map the input coordinate to a higher dimensional space through \(positional\)\(encoding\) before feeding them to the network. Successful implementations of \(positional\)\(encoding\) in seismic applications can be found in the works of [19] and [20]. In our implementation, we also leverage positional encoding to ensure that the network is capable to reconstruct multi-scale, oscillating signals like those encountered in seismic data.

To support our claims, we present numerical examples focused on two of the most challenging tasks in seismic interpolation: namely, interpolation of regularly subsampled data (beyond aliasing) and data with large-gaps. We first evaluate the performance of the proposed approach on two simple synthetic data examples, comparing the results obtained by the PINNs approach with a simple plane-wave regularized least-square inversion (PWLS) and against a previous version of our framework called PWD-PINN ([21]). We then consider a field data example where given the challenge of the higher frequency content and higher complexity of the recorded signals, we decided to leverage on positional encoding to overcome the low frequency bias of these types of architectures ([18]). The proposed framework, that we refer to as PINNslope, not only allows to interpolate seismic data, but it can also be used to estimate slopes from fully sampled data of quality comparable to those of the PWD filters. In the context of sparsely sampled data, the procedure is advantageous because we can directly estimate the slopes during the interpolation of the recorded wavefield (i.e. using the aliased data), without the need of low-pass filtering the data to obtain an alias-free version on which to perform a reliable slope estimation as in our previous work. The inversion result of the slope network turn outs to be smooth, but a more accurate version than the one estimated from the low frequency data by means of PWD filters. Thanks to that the physical constrain the network is very efficient in steering the interpolation process towards an accurate physically driven solution.

To summarize, our main contributions comprise:

1. A novel machine learning framework for slope assisted seismic data interpolation.
2. An innovative procedure for local slope attribute estimation by the mean of physics informed neural networks.
3. Successful application of the framework on field data.

The paper is organized as follows. First, we present the theoretical background of our methodology. We then describe the network architecture, focusing on some of the key implementation details needed to achieve a stable training process. Finally, our method is applied to a range of synthetic and field data. Its reconstruction capabilities are compared to those of a classic least-squares inversion regularized by the discretised local plane-wave equation (PWLS), as well as to our previous published PWD-PINN approach which uses pre-computed local slopes.
Theory

### Problem statement

The objective of this paper is to formulate the problem of seismic data interpolation within the framework of physics informed neural networks.

To begin with let's define the basic mathematical model used to obtain a decimated version of the original seismic data: a restriction operator **R**, is defined such that it samples the columns of the data matrix \(\textbf{u}=[\textbf{u}_{1}^{T},\textbf{u}_{2}^{T},...,\textbf{u}_{N}^{T}]\) at desired locations (where \(N\) is the number of total traces in the dataset), removing the missing traces that we wish to interpolate. In matrix notation the operation of subsampling a gather of traces can be written as:

\[\textbf{d}=\textbf{Ru} \tag{1}\]

where **d** is the subsampled data with either missing traces at a regular interval to simulate spatial aliasing, or missing a large number of consecutive traces to simulate a gap in the acquisition geometry.

### Slope estimation with plane-wave destructors

The physical model used to express seismic data as local plane-waves is represented by the local plane-wave differential equation:

\[\frac{\partial u(t,x)}{\partial x}+\sigma(t,x)\frac{\partial u(t,x)}{\partial t }=r(t,x)\approx 0, \tag{2}\]

where \(u(t,x)\) is the pressure wavefield, \(r(t,x)\) is the PDE residual, and the parameter \(\sigma(t,x)\) is the local slope (or wavenumber) with units equal to the inverse of the velocity of propagation. An analytical expression exists for the solution of equation 2 in case of a constant slope, which is simply represented by a plane-wave

\[u(t,x)=f(t-\sigma x), \tag{3}\]

where \(f(t)\) is an arbitrary waveform at \(x=0\). We can see that the left hand side of equation 2 decreases as the observation \(u(t,x)\) matches the wave displacement \(u(t-\sigma x)\) ([22]).

In our work, we are interested in computing a slope varying both in time and space, but no analytical solution exists for such a case. Hence, Claerbout (1992) casts the dip estimation as a linear least-squares problem, through an operation named plane-wave destruction. In this approach, the curvature of the events is linearly approximated by computing the slope in a small window of the entire data. The slope is estimated through equation 2, minimizing the quadratic residual:

\[Q(\sigma)=(\textbf{u}_{x}+\sigma\textbf{u}_{t})\cdot(\textbf{u}_{x}+\sigma \textbf{u}_{t}), \tag{4}\]

where \(\textbf{u}_{x}\) and \(\textbf{u}_{t}\) are respectively defined as the spatial and temporal derivatives of the wavefield **u**. Setting the derivative of \(Q(\sigma)\) to zero, we can find its minimum as:

\[\sigma=-\frac{\textbf{u}_{x}\cdot\textbf{u}_{t}}{\textbf{u}_{t}\cdot\textbf{ u}_{t}}. \tag{5}\]

Fomel (2002), on the other hand, frames the slope estimation as a non-linear least-squares problem computing the slope globally on the entire data by the means of a plane-wave destruction filters. Given a wavefield **u** as a gather of seismic traces \(\textbf{u}=[\textbf{u}_{1},\textbf{u}_{2},...,\textbf{u}_{N}]^{T}\), the destruction operator predicts each trace from the previous one by shifting the observed trace along the dominant local slopes of the seismic data and subtracts the prediction from the original one. The phase-shift operation on the traces is approximated by an all-pass digital filter (or prediction filter in 2D). The filter coefficients are determined by fitting the filter frequency response (at low frequencies) to the response of the phase-shift operator. In this implementation, the slope (\(\sigma\)) enters in the filters coefficients in a non linear way. To characterize and describe the entire gather, the prediction of several plane waves is needed (and not only one). This is achieved by cascading various filters of the above mentioned form. The filter is applied to the data **u** as a convolutional operator \(\textbf{D}(\sigma)\). In matrix notation, the slope estimation problem can be written as

\[\textbf{D}(\sigma)\textbf{u}=\textbf{r}, \tag{6}\]
where **r** is the residual. This non-linear least-squares problem is solved via Gauss-Newton iterations, which implies solving

\[\textbf{D}^{\prime}(\sigma_{0})\Delta\sigma\textbf{u}+\textbf{D}(\sigma_{0}) \textbf{u}=\textbf{r}, \tag{7}\]

where \(\textbf{D}^{\prime}(\sigma_{0})\) is the derivative of the filter coefficients \(\textbf{D}(\sigma)\) with respect to \(\sigma\). The minimization problem is solved for the update \(\Delta\sigma\), which is repeatedly added (at every iteration) for an initial guess \(\sigma_{0}\). The problem can be regularized by adding an appropriate penalty term that avoids oscillatory solutions of the slope attribute.

### Plane-wave regularized least-squares interpolation

In this section, we describe a conventional approach to take into account pre-computed slopes whilst interpolating seismic data (i.e., restoring missing traces). This method will be later used as a benchmark for our PINNs approach. The inverse problem is cast as follows: finding the shot-gather **u** (i.e., the full gather of traces) that minimizes the Euclidean distance between the subsampled data **d** and the estimated subsampled data **Ru**, whilst at the same time satisfying the plane-wave differential equation with pre-computed slopes. The objective function is formally defined as

\[f(\textbf{u})=\|\textbf{d}-\textbf{Ru}\|_{2}^{2}+\epsilon_{r}\|\textbf{u}_{x }+\boldsymbol{\Sigma}\textbf{u}_{t}\|_{2}^{2} \tag{8}\]

where \(\textbf{u}_{x}\) and \(\textbf{u}_{t}\) are respectively defined as the spatial and temporal derivatives of the data **u**, \(\boldsymbol{\Sigma}\) is a diagonal matrix that applies element-wise multiplication of the pre-computed local slope and \(\epsilon_{r}\) is a weight to control the contribution of the PDE in the solution. The data term of the objective function aims at accurately reproducing the available traces from the estimated full shot-gather **u** subsampled by the restriction operator **R**. This means that all the interpolation operation between the traces is performed by the local plane-wave regularization term. In other words, the regularization term has the function of filling the gaps between the subsampled traces, spraying the information available from two neighboring traces along the curvature of the provided local slope field.

### Physics Informed Neural Networks

In this section, we aim to show that starting from a knowledge of a slope field, estimated via PWD filters or any other algorithm, the problem of seismic interpolation can be formulated within the PINNs framework.

PINNs have been designed to blend the universal function approximator capabilities of neural networks ([23]) with a physical constraint given by a PDE, which describes the physical system under study. In our specific case, the PDE that we seek to satisfy is the local plane-wave differential equation (equation 2).

A neural network \(\phi_{\theta}(t,x)\) is designed to approximate the function \(u(t,x)\), where \(\theta\) refers to the weights (and biases) to be optimized and the pair (\(t,x\)) represents the input to the network. The network predicts the recorded wavefield \(u(t,x)\) at the corresponding location in the time-space domain of interest. A remarkable convenience of PINNs is that in contrast to traditional numerical methods, they do not require a discretization of the computational domain. The partial derivatives of the underlying PDEs are computed by the means of automatic differentiation (AD), which is a general and efficient way to compute derivatives based on the chain rule. AD is usually implemented in neural networks training to compute the derivatives of the loss function with respect to the parameters of the network. However, AD can be more broadly applied to every computational program that performs simple arithmetic operations and calculates elementary functions (linear transformations and non-linear activation functions in the case of neural networks) by keeping track of the operations dependencies via a computational graph and successively computing their derivatives using the chain rule. The PINN framework is trained in an unsupervised manner, using a loss function which includes both the local plane differential equation and a set of \(N_{t}\) (number of traces in the subsampled gather) boundary conditions corresponding to the available traces

\[\mathcal{L}=\frac{1}{N_{u}}\sum_{i=0}^{N_{u}}\left(\frac{\partial\phi_{\theta }(t_{i},x_{i})}{\partial x}+\sigma(t_{i},x_{i})\frac{\partial\phi_{\theta}(t_ {i},x_{i})}{\partial t}\right)^{2}+\lambda\left(\frac{1}{N_{t}}\sum_{j=0}^{N_ {t}}|u(t_{j},x_{j})-\phi_{\theta}(t_{j},x_{j})|\right), \tag{9}\]

where \((t_{i},x_{i})\) are points randomly sampled from the input space with \(N_{u}\) as number of total grid points, \(u(t_{j},x_{j})\) is the known solution regions at points indexed by \(j\) (available traces), and \(\lambda\) is a scalar weight for the second term.

In this first approach named PWD-PINN (1), the slope \(\sigma(t_{i},x_{i})\) is pre-computed and it is estimated by means of PWD filters, outside the training process of the network. The slope array remains fixed during training and it is not updated.
### Simultaneous data interpolation and slope estimation

In this section, we introduce the slope estimation framework using physics informed neural networks, named PINNslope. We propose to estimate the local slopes while at the same time interpolating the aliased data (or any other type of interpolation task). Specifically, we simultaneously train two neural networks to predict the data and the local slopes that satisfy the plane-wave PDE. This approach bears similarity with previous works by [12] and [24], in the context of traveltime tomography.

As shown in the diagram of fig.2, both networks have fully-connected architectures and utilize \(Tanh(\cdot)\) activation functions; moreover, a positional encoding is added to the wavefield network to tackle multi-scale signals (i.e., signal exhibiting both low and high frequency components) such as our seismic traces. The two networks also differ in the number and size of the layers: the wavefield network aims at reconstructing the shot gather data and requires a much larger number of degrees of freedom to fit the complexity of seismic signals (i.e., the seismic traces \(u(t_{j},x_{j})\)); rather for the slope, as we aim to obtain a smooth solution, this can be achieved by using more compact architecture. After computing the loss function, both networks are simultaneously updated. Two separate ADAM optimizers are utilized to allow two different learning rate values if necessary.

### Positional Encoding

During the numerical experiments both frameworks struggled in fitting signals with high frequency content. In our previous experiments, the low frequency bias of neural networks ([18]) was addressed with frequency upscaling by the

Figure 1: A diagram explaining the PWD-PINN algorithm. The network is trained while maintaining the slope array fixed during training.

Figure 2: Double network scheme for joint wavefield and slope estimation in the same training procedure.
mean of _neuron splitting_ ([25]) and with locally adaptive activation function ([26]). Here, the low frequency bias of multi-layer perceptrons (MLPs) is tackled by including _positional encoding_ of spatial coordinates ([16]).

Differently from the classical Transformer approach to positional encoding, where it is used to track the token positions, in our application we use it to map the input coordinate grid into a higher dimensional space, which allows for a better fit of high frequency signals. The approach implemented in this work resembles the one previously presented in [17] and referred to as _Fourier feature mapping_, where the authors utilized a linear sampling in the Fourier space that enables for a large amount of frequency components in the low-frequency regions. This modification is fundamental also for our implementation as other forms of encodings were introducing noise into the reconstruction, as the higher frequencies present in the encodings were able to fit the noise in the traces. Additionally, [20] proposed an anisotropic version of positional encoding justified by the idea that seismic data components present different features and should not be equally encoded. The formulation utilized here can be summarized as follows:

\[\gamma_{X}(x)=[cos(k_{0}x),sin(k_{0}x),...,cos(k_{X-1}x),sin(k_{X-1}x)], \tag{10}\]

where \(k_{x}=\frac{\pi x}{2}\) with \(x=0,...,X-1\), \(\in\mathbb{N}\), is a simple linear sampling, and \(X\) represents the number of encoded frequencies in the \(x\) coordinate. The \(t\) coordinate is encoded in the same way with the number of frequencies equal to \(T\)\(\in\mathbb{N}\), and the encoded coordinates are subsequently concatenated together as

\[\Gamma_{X,T}(x,t)=[\gamma_{X}\left(x\right),\gamma_{T}(t)]. \tag{11}\]

The positional encoding operation has been embedded as a network layer inside the architecture of the data network, and the number of frequencies corresponding to each coordinate is decided through trial and errors.

## 3 Numerical experiments

In this section, the proposed methodology is tested on synthetic and field data. For both the PWDPINN and PINNslope approaches, a feed-forward neural network architecture with 4 layers and a \(Tanh(\cdot)\) activation function was utilized. We set the number of neurons to be the same for all layers in all our experiments but this number can vary in different experiments as we specify it in the subsections. In both frameworks, the networks are trained in an unsupervised manner, passing as input an ensemble of \((x,t)\) points. The ensemble is passed to networks in batches of 1000 randomly sampled points. For every batch, the ensemble of the collocation points is concatenated to an array containing half of the points \((x_{j},t_{j})\) associated with the available traces to be fitted. All networks in every experiments are trained using ADAM optimizer, with learning rate fixed at \(10^{-3}\). These parameters are chosen based on some initial tests and kept fixed throughout the study.

### Synthetic data examples

#### 3.1.1 Local slope estimation

In this first example, we estimate the slope with the PWD algorithm and with the PINNslope framework, to compare their performance. The synthetic seismic image (Sigmoid model, [1]) is assumed to be fully sampled and all the traces have been utilized in the training process.

Figure 3: Slope comparison between PWD algorithm. a) Seismic image, b) local slope estimate from Plane-wave destruction filters [2], c) local slope estimate obtained with the PINNslope approach, d) difference between the PWD and PINNslope estimated slopes.
As shown in Fig.3, the PINNslope framework can accurately estimate the local slope of complex subsurface geometries, and it results in a slightly smoother version with less artefacts near the major fault compared to the local slope estimated via the PWD algorithm.

#### 3.1.2 Interpolation beyond aliasing with local slope estimation

The goal of this second example is to reduce the spatial aliasing present in the recorded data by interpolating the missing traces. The synthetic data in fig.4a have a trace spacing of 10 meters and have been subsampled by a factor of 5 through the operator **R**, to obtain the aliased version in fig.4b.

It is not possible to apply the PWD filter directly to estimate the slope from the subsampled data, since it can lead to erroneous estimates by picking the aliased dips instead of the true ones of the fully sampled data [2].

To avoid this issue, the following pre-processing steps have been performed:

1. Apply \(f-k\) filter to the spectrum of the aliased data.
2. Inverse transform the filtered spectrum to get a low frequency alias-free version of the data.
3. Apply the PWD filters algorithm to the low frequency data and estimate the slope.
4. Utilize the PWD estimated slope from low frequency data inside the PINN loss function.

The network capacity corresponds to 4 layers with 512 neurons each, with the number of encoding frequencies set to \(X=8\) and \(T=32\) for the \(x\) and \(t\) coordinates respectively. The network is trained using the loss function in eq.9 with the parameter \(\lambda\) set to 1000 and \(\sigma\) correspond to the PWD estimated slope displayed in fig. 4c.

Fig.5 compares the results obtained with the different approaches. The output of the regularized least-squares inversion in fig.5a demonstrates the importance of the plane-wave penalty term, which helps in filling the gap between the

Figure 4: a) Original seismic data, b) seismic data with missing traces, c) low frequency data from which the PWD slope has been computed, d) \(f-k\) spectrum of the data, e) \(f-k\) spectrum of the subsampled data, f) PWD estimated slope from the low frequency data in fig.4c.
available traces following the correct overall geometry of the arrivals. Unfortunately, as soon as the reflections start bending their resolution decrease, worsening towards the far-offset. In this interpolation attempt the sharp and definite seismic response that characterize this simple synthetic data is slightly spread in a fuzzy pattern, a sign that the algorithm cannot properly restore the energy in the correct position. The difference with the original data in fig.5d shows the amount of energy lost, as well as some artifacts. The achieved result is almost perfect where the arrivals are generally linear. Anyhow, LS inversion is almost instantaneous compared to the neural networks approach. The result displayed in fig.5b requires a runtime of approximately 26 minutes for 2000 epochs as shown in the plot of the loss curves in fig.6a.

As in the previous result, the quality of the interpolation of the PWD-PINN algorithm decreases in the far-offset, although only for the first few reflections. This is a limitation of the algorithms that leverage on the PWD estimated slope, which is inaccurate at the far-offsets where the events are steeper; the slope estimated via PWD inherently contains errors because of the procedure through which it has been computed, but even more it has been estimated from a low frequency version of the original data. Despite the poor interpolation of the above mentioned arrivals, all the others look adequately restored. Most of the energy is in the correct position as we can see from its spectrum in fig.5h.

Figure 5: a) Plane-wave regularized least-squares inversion interpolation result, b) PWD-PINN interpolation result, c) PINNslope interpolation result, d) difference between fig.4a and fig.5a, e) difference between fig.4a and fig.5b, f) difference between fig.4a and fig.5c, g) \(f-k\) spectrum of fig.5a, h) \(f-k\) spectrum of fig.5b, i) \(f-k\) spectrum of fig.5c.
Anyhow, in this result the resolution is lower, in fact the traces interpolated at the far-offsets include gaps. Moreover, in the first two events the amplitude is not properly reproduced.

The best reconstruction is clearly given by the PINNslope framework. The architecture of the network has the same capacity then the one of the PWD-PINN algorithm and the loss function in fig.6b shows that it has been trained for the same amount of epochs as PWD-PINN. The key difference in the result is made by the second smaller network that approximates the local slope function. The slope estimate is carried on simultaneously with the interpolation performed by the bigger data network on the original shot gather (Fig. 4a), no filtering required. This simultaneous updating process of data and slope allows for larger search space to speed up convergence to the accurate data. As can see from fig.7, the PINN estimated slope closely matches the accurate PWD slope computed from the data of figure 4a (ignoring the right-upper part where there are no arrivals and the two approaches clearly extrapolate the values).

### Field data examples

#### 3.2.1 Interpolation beyond aliasing with local slope estimation

The numerical examples below are performed on a field dataset from the Gulf of Mexico. Here, only the results for the PINNslope approach and PWD-PINN will be compared as we want to focus on the PINN approaches. The trace spacing in the original shot-gather of fig.8a is 26.7 meters and we subsampled it by a factor of 5, increasing the spacing between the traces to 133.5 meters (fig.8b).

Figure 6: a) Loss of the PWD-PINN training. b) Loss resulting from the PINNslope training. In the small box inside the two plots it is shown the contribution of each term of the loss function in log-scale, that is: data-term (in green), physical-term (in blue) and total loss (in orange).

Figure 7: a) Accurate slope estimated with PWD algorithm from the data shown in 4a. b)Slope estimate through the PINNslope framework while simultaneously interpolating the data in fig.5c, c) difference between fig.7a and fig.7b.
As mentioned earlier, computing directly the slope from the subsampled gather is not feasible and some pre-processing steps are required. Filtering out the aliased part of the field data in fig.8e is far more challenging then for the on synthetic ones; the part of the signal that is not aliased is very small and does not contain significant energy. The retrieved low frequency data (fig.8c) are fed into the PWD algorithm and, due to the low frequency nature of the data, the resulting local slope is a low resolution rough estimate of the slope of the high-resolution data (fig.8f). The network has an architecture of 4 layers and 512 neurons in each layer. It is equal to the one used for synthetic data. That is because even if the traces are way more complex in the field data, the synthetic traces have an amplitude that does not decrease as much in time. From our initial tests the network requires the same capacity to easily fit the strong oscillations of the synthetic signals. The slope network (as in the synthetic case) has 2 layers with 2 neurons to estimate a smooth version of the slope field.

The PDW-PINN algorithm (fig.9a) does not achieve a good result. Only the interpolation in the near-offset could be consider reasonable. In the far-offsets the reconstruction is worse as we already observed in the synthetic data. In this part of the dataset, the PWD estimated slope is prone to errors and it does not allow a good interpolation. The loss curves (fig.10a) show that we had to increase the \(\lambda\) parameter to very high values (\(\lambda=10000\)) to make the network properly fit the traces. The network struggles to accurately fit the traces if the accuracy of the local slope is poor, as this will negatively effect the PDE term of the loss function. As a result, increasing the weight on the data fitting term rendered the PWD-PINN algorithm to be a data fitting algorithm.

In contrast, the PINNslope approach admits a good performance. It reproduces the original dataset; most of the energy has been restored and the aliasing has been suppressed. The extra degrees of freedom provided by the small slope network helped the convergence.

Figure 8: a) The original seismic data, b) subsampled seismic data, c) low frequency data from which the PWD slope has been computed, d) \(f-k\) spectrum of the original data, e) \(f-k\) spectrum of the subsampled data, f) PWD estimated slope from the low frequency data in fig.8c.
In fig.11, the PINN estimated slope is compared to the PWD slope computed from the full data (fig.8a). The PINN slope is smoother than the PWD one and again admits generally lower values, probably because it has been estimate on less dense data. However, the overall trend of the PINNslope slope is correct and its smoothness serves to its purpose in the plane-wave regularization term. If compared to the PWD slope in fig.8f, which is the realistically achievable slope when we try to solve an interpolation problem of this kind, the PINN slope is a way better and more precise estimate.

The residuals shown in fig.9f and fig.9e, differently from the synthetic data, are partially due to the field data containing secondary events with conflicting dips that cannot be recovered by our method. We note that this is a general weakness of interpolation methods relying on the plane-wave PDE.

#### 3.2.2 Performance assessment

In this section, the PINNslope framework is tasked to a harder interpolation where fewer traces are available. The aim is to assess its performance on the current dataset and evaluate its limits. Moreover, as subsequent shot gathers in the dataset will have only minor changes between them, we test the converge behaviour of the pre-trained network when applied to the next gathers is the data. We first apply PINNslope on the shot-gather subsampled by various factors: 6, 7 and 8 (respectively, 160.2 meters, 186.9 meters, 213.6 meters intervals between each traces). Finally, we also test the ability of our framework to interpolate a dataset with a large gap of traces (i.e. 15 traces, for a total of 400.5

Figure 9: a) PWD-PINN interpolation result, b) Difference between the original full data (fig.8a) and PWD-PINN result, c) \(f-k\) spectrum of PWD-PINN result, d) PINNslope interpolation result, e) difference between the original full data (fig.8a) and PINNslope result, f) \(f-k\) spectrum of PINNslope result.
meters gap) placed in the middle of the gather. This a very hard task for the network, which has to rely solely on the information obtainable from the left and right side of the gather, as in the gap region, there is no knowledge of the shape of the arrivals nor no that of their slope field. So far, we are not aware of any interpolation algorithm that can solve this category of interpolation tasks in an automatic and physically driven manner.

The performance results are shown in fig.12. For this shot gather, a higher subsampling of 6 and 7 (respectively fig.12a and fig.12b) does not impact the network performance. The signal to noise ratio is maintained almost constant from the result described in the previous section and the arrivals are perfectly interpolated. The framework start to face some challenge when the gather is subsampled by a factor of 8 (213.6 meters interval between the traces). In the near-offset the interpolation is still accurate but as the dip of the arrival starts to increase, PINNslope is unable to retrieve the correct slopes and struggles in interpolating the arrivals. We consider this sampling to be the threshold limit for the framework for this gather (and the freqeuncy range involved). The presented result are obtained with an increasing the number of epochs as the subsampling increases (1500 epochs for a subsampling of 6, and 2500 for a subsampling of 8). In the case of the higher subsampling factor, it was additionally needed a \(\lambda\) value of 10000 instead of the usual 100.

By enlarging fig.12 at early times (fig.13) we can see the main difference between the original data and the PINNslope interpolation. In fact, most of the events are well reproduced, as what is missing is the energy corresponding to the often weaker events with conflicting; as such the interpolated data looks cleaner than the original one. As mentioned before most algorithms that rely on slope estimation to perform their tasks cannot leverage on the energy of the second

Figure 11: a) The slope estimated with PWD algorithm from the data shown in 8a. b) The slope estimate through the PINNslope framework while simultaneously interpolating the data in fig.9e, c) the difference between fig.11a and fig.11b.

Figure 10: a) Loss curves of the PWD-PINN training on real data. b) Loss curves of the PINNslope training on real data. In the small box inside the two plots we show the contribution of each term of the loss function in log-scale, that is: data-term (in green), physical-term (in blue) and total loss (in orange).
order dips for its reconstruction unless some additional slopes are included in the process.

The result of the gap interpolation (fig.12d and fig.13c) is remarkable. The PINNslope framework is able to extrapolate the main reflection from the left and right and connect them together. Of course we are aware of the high errors in the middle of the gap compared it to the same part in the original data. The algorithm cannot completely restore the missing part, as mentioned above it lacks information on the conflicting dips. Is it worth mentioning that this experiment has been done to show the interpolation capabilities of neural networks and especially of PINNs. Fig.14 shows the convergence capabilities of the PINNslope network (pre-trained on the gather of fig.9a with a subsampling factor of 5) on the subsequent gathers of the dataset (that have been subsampled by a factor of 5).

Specifically, the gathers of fig.14e and fig.14f are consecutive to the one used as benchmark until now. Instead, the ones in fig.14g and fig.14h are way more distant. The idea is to use only one pre-training step (on the gather of fig.9a) and a small fine-tuning step to make the network fit the subsequent gather. As we can see, the number of epochs needed for fine-tuning the network depend on how different the gathers are from the one used for the pre-training. Although, all the gathers has been well reproduced with a small number of epochs.

Figure 12: The performance of PINNslope with high subsampling rates and a particular case where an entire part of the gather is missing: a) the result obtained starting from data subsampled by 6, b) result obtained starting from data subsampled by 7, c) result obtained starting from data subsampled by 8, d) result obtained starting from data that contains a large gap of 15 traces (white hatched lines delimit the gap area).

Figure 13: Zoom in of some of the results in fig.12. a) Zoom in of the original data, b) Zoom in of the result obtained starting from a data subsampled by a factor of 7, c) Zoom in of the result starting from a data that contains a large gap of 15 traces (white hatched lines delimit the gap area).
## 4 Discussion

Out of all the implementations we tested, PINNslope achieves the best performance. It can interpolate the missing data eliminating the alias present in the data whether they are related to coarse recordings or the presence of obstacles during the acquisition. Parameterising the plane-wave PDE in the loss function with a small neural network does not significantly affect the runtime with respect to PWD-PINN. For PWD-PINN and PWLS inversion a number of pre-processing steps should be performed before being able to compute the local slope to use in the PDE regularization term. This is also a time consuming part of these approaches. Although, the PWLS inversion is almost instantaneous compared to the PINNs.

All of the three methodologies present similar drawbacks due to the of the plane-wave approximation. In fact, the estimated slopes with the PWD filters (or any other algorithm that can estimate slopes) and PINNslope are computed with respect to the main events and they are not able to retrieve information for the conflicting slopes (slopes that have opposite or different direction) in the data. The PINNslope framework cannot reproduce the events that have a reverse slope with respect to the main trend of the arrivals. Moreover, the estimate of the slope on a data point where two events are crossing each other will always induce errors in the data fit. Currently, it is not possible to compute the value of two slopes on one single point of the dataset unless we include two slopes in the framework, which will be investigated in the future.

PINNslope has the potential to fill a large gap of traces via a physically driven approach not achievable by any other type of algorithm, furthermore estimating the full slope field. This example illustrates the interpolation capabilities of this type of implementation. A substantial help to the network fitting ability is derived from the _positional encoding_ layer. In the past experiments and implementations various methodologies have been tested to make the network able to fit complex signals as field seismic traces. _locally-adaptive activation functions_ ([26]) has been implemented along with a _Sin()_ activation function ([21]), to allow the network to be more expressive without the need to extend its capacity, but still it was not enough. In addition, a technique named frequency upscaling by the mean of neuron splitting ([25]) have been tested achieving good results in fitting complex high frequency signals. One of its drawbacks is that it was injecting low level noise in the reconstructed result and required to train various time the network to adequately fit the entire frequency content of the dataset. On the other hand, _positional encoding_ solved this issue allowing the network to

Figure 14: The results of fine-tuning PINNslope on different gathers from the one used to pre-train the network (fig.9a). a) interpolated gather subsequent to the one in fig.9a, b) interpolated gather subsequent to the one in fig.14b, c) interpolated gather at position 100 in the dataset, d) interpolated gather at position 150 in the dataset, e) original gather at position 2, f) original gather at position 3, g) original gather at position 100, h) original gather at position 150.
accurately reproduce the field data and achieve a faster convergence, without the need of training several times on the same dataset (as frequency upscaling with neuron splitting).

## 5 Conclusions

We introduced a novel PINN framework in the seismic signal processing field for simultaneous seismic data interpolation and local slope estimation. It is discernible also as an innovative procedure for local slope attribute estimation of complex subsurface images. The obtained results are compared on synthetics and field datasets against PWD-PINN, which showed good results previously in the synthetic data experiments, and the plane-wave regularized least-square inversion (for the synthetic example) to better examine the PINNs performances against an approach that does not rely on neural networks to optimize the same objective function.

We found that introducing a second network to estimate the local slope attribute while at the same time interpolating the aliased data achieves better results in terms of signal to noise ratio, while also improving the overall network convergence. The positional encoding layer was a fundamental addition to the architecture and helped overcome previous difficulties such as high frequency fitting and noise introduction during the interpolation process. The PINN estimated slopes look accurate and consistent with the interpolated data, their accuracy is comparable to the one obtainable with the plane-wave destruction filters estimate.

## Acknowledgments

This publication is based on work supported by the King Abdullah University of Science and Technology (KAUST). The authors thank the DeepWave sponsors for supporting this research.
