Deep Convolutional Neural Network and Sparse Least Squares Migration

TL;DR

This paper links deep CNNs to seismic migration, proposing a neural network version of sparse least squares migration that reduces computation and enhances noise removal for high-resolution imaging.

Contribution

It introduces neural network sparse least squares migration (NNLSM), connecting CNN filters to seismic Green's functions, and offers a computationally efficient, noise-robust imaging method.

Findings

NNLSM reduces computational cost compared to standard LSM.

NNLSM effectively denoises seismic migration images.

Quasi-reflectivity images serve as superresolution attributes.

Abstract

We recast the forward pass of a multilayered convolutional neural network (CNN) as the solution to the problem of sparse least squares migration (LSM). The CNN filters and feature maps are shown to be analogous, but not equivalent, to the migration Green's functions and the quasi-reflectivity distribution, respectively. This provides a physical interpretation of the filters and feature maps in deep CNN in terms of the operators for seismic imaging. Motivated by the connection between sparse LSM and CNN, we propose the neural network version of sparse LSM. Unlike the standard LSM method that finds the optimal reflectivity image, neural network LSM (NNLSM) finds both the optimal quasi-reflectivity image and the quasi-migration Green's functions. These quasi-migration-Green's functions are also denoted as the convolutional filters in a CNN and are similar to migration Green's functions. The advantage of NNLSM over standard LSM is that its computational cost is significantly less and it can be used for denoising coherent and incoherent noise in migration images. Its disadvantage is that the NNLSM quasi-reflectivity image is only an approximation to the actual reflectivity distribution. However, the quasi-reflectivity image can be used as a superresolution attribute image for high-resolution delineation of geologic bodies.

Equations

$$Given:$$

$$Find:$$

$${\mathbf{m}}^{\ast }$$

$${\mathbf{m}}^{(k+1)}$$

$$ϵ$$

$${\mathbf{m}}^{\ast }$$

$$m_i^{(k+1)}$$

$$\displaystyle m_{i}^{(k+1)}={\mbox{soft}}({\Big{[}}{\bf m}^{(k)}-\frac{1}{\alpha}{\bm{\Gamma}}^{T}({\bm{\Gamma}}{\bf m}^{(k)}-{\bf m}_{mig}){\Big{]}}_{i},\frac{\lambda}{\alpha}),$$

$${\mathbf{m}}^{(1)}=\text{mbox}soft({\mathbf{\Gamma }}^T{\mathbf{m}}^{mig},\lambda ).$$

$${\mathbf{m}}^{(1)}={\sigma }_{\lambda }({\mathbf{\Gamma }}^T{\mathbf{m}}^{mig}).$$

$${\mathbf{m}}^{mig}=\sum _i^{k_0}{\mathbf{\gamma }}_i\ast {\mathbf{m}}_{1i}^{\prime }.$$

$${\mathbf{m}}_i,{\Gamma }_i\text{mbox}suchthat$$

$${\mathbf{m}}_1^{\ast }$$

$${\mathbf{m}}_2^{\ast }$$

$${\mathbf{m}}_N^{\ast }$$

$${\mathbf{\Gamma }}^{(i)}={\mathbf{\Gamma }}_1{\mathbf{\Gamma }}_2\dots {\mathbf{\Gamma }}_i,$$

$${\mathbf{m}}^{mig}={\mathbf{\Gamma }}^{(i)}{\mathbf{m}}_i.$$

$$\text{mbox}for\mathbf{x}\in D_{model}:m(\mathbf{x})=$$

$$\eta {\int }_{D_{model}}d\mathbf{z}\stackrel{\Gamma (\mathbf{x}|\mathbf{z})}{\overbrace{{\int }_{\mathbf{y}\in D_{data}}d\mathbf{y}{\omega }^{4}G(\mathbf{x}|\mathbf{y}{)}^{2\ast }G(\mathbf{y}|\mathbf{z}{)}^2}}m(\mathbf{z}),$$

$$\text{mbox}for\mathbf{x},\mathbf{z}\in D_{model}:\Gamma (\mathbf{x}|\mathbf{z})=$$

$$\eta {\int }_{\mathbf{y}\in D_{data}}d\mathbf{y}{\omega }^{4}G(\mathbf{x}|\mathbf{y}{)}^{2\ast }G(\mathbf{y}|\mathbf{z}{)}^2.$$

$$[\mathbf{\Gamma }\mathbf{m}{]}_i$$

$$ϵ$$

$$\frac{\partial ϵ}{\partial x_i}$$

$$\frac{\partial ||\mathbf{x}|{|}_1}{\partial x_i}=1for{x}_i\ge 0;$$

$$x_{i}=\mbox{soft}(z_{i},\lambda)={\left\{\begin{array}[]{cl}z_{i}-\lambda&if~{}~{}z_{i}~{}~{}\geq~{}~{}\lambda\\ 0&if~{}~{}|z_{i}|<~{}~{}\lambda\\ z_{i}+\lambda&if~{}~{}z_{i}~{}~{}<-\lambda\end{array}\right\}}.$$

$${\mathbf{x}}^{\ast }$$

$$x_i^{(k+1)}$$

$$\displaystyle(\tilde{\bf m}^{*},\tilde{\bm{\Gamma}}^{*})=\arg\min_{\tilde{\bf m},~{}\tilde{\bm{\Gamma}}}{\Big{[}}\frac{1}{2}||\tilde{\bm{\Gamma}}\tilde{\bf m}-{\bf m}^{mig}||_{2}^{2}+\lambda S(\tilde{\bf m}){\Big{]}},$$

$$\displaystyle\epsilon=\frac{1}{2}\sum_{i}{\Big{[}}\sum_{j}\tilde{\Gamma}({\bf x}_{i}|{\bf z}_{j})\tilde{m}_{j}-m_{i}^{mig}{\Big{]}}^{2}+\lambda S(\tilde{\bf m}).$$

$$\frac{\partial ϵ}{\partial \tilde{\Gamma }({\mathbf{x}}_{i^{\prime }}|{\mathbf{z}}_{j^{\prime }})}$$

Full text

1 Introduction

Deep convolutional neural networks (CNNs) have been recently used for solving geophysical problems, such as seismic first-arrival picking (Lu and Feng,, 2018; Yuan et al.,, 2018; Hu et al.,, 2019), seismic interpretation (Wu et al., 2019a, ; Wu et al., 2019b, ; Shi et al.,, 2019), and seismic imaging and inversion (Xu et al.,, 2019; Sun et al.,, 2020; Kaur et al.,, 2020). Interpretation examples include salt classification (Waldeland et al.,, 2018; Shi et al.,, 2019), fault detection (Huang et al.,, 2017; Xiong et al.,, 2018; Wu et al., 2019b, ; Zheng et al.,, 2019), reservoir characterization (Karimpouli et al.,, 2010; Cao and Roy,, 2017; Zhu et al.,, 2017), and seismic lithofacies classification (Ross and Cole,, 2017; Liu et al.,, 2019) with semi-supervised learning (Di et al.,, 2020). Other applications include the use of neural networks for well-log interpolation (Saggaf and Nebrija,, 2003; Salehi et al.,, 2017; Pham et al.,, 2020), seismic data interpolation (Mandelli et al.,, 2018, 2019; Wang et al.,, 2020), velocity model building (Araya-Polo et al.,, 2018; Richardson,, 2018), well-log ties (Bader et al.,, 2019), synthetic well-data generation (Rolon et al.,, 2009), autoencoders for unsupervised facies analysis (Qian et al.,, 2018), and supervised horizon tracking (Peters et al., 2019a, ; Peters et al., 2019b, ). Recently, an unsupervised autoencoder method with regularization was developed by Shi et al., (2020) to track target horizons.

There are many types of neural network or machine learning methods, selections ranging from generative adversarial networks for seismic interpolation (Siahkoohi et al.,, 2018), residual networks for traces missing at regular intervals (Wang et al.,, 2019), Monte Carlo and support vector regression (Jia et al.,, 2018) for data interpolation, autoencoders (Wang et al.,, 2020; Shi et al.,, 2020) for target horizons tracking, and recurrent neural networks for well-log interpolation (Pham et al.,, 2020). The design of geophysical CNN architectures has largely been based on empirical evidence from computer vision research, insights from the principles of artificial intelligence and heuristic experimentation. Heuristic experimentation is most often used to decide the parameter design111The design of a CNN architecture selects the number of CNN layers, the number of filters/layer, the size of the filters, the type of activation functions, the number of skip layers and whether a layer acts as a decoding or encoding operation. for the CNN architecture, which has both merits and liabilities. The merit is that trial-and-error with different architecture parameters is likely to give excellent results for a particular data set, but it might not be the best one for a completely different data set. This shortcoming in using empirical tests for parameter selection largely results from the absence of a rigorous mathematical foundation (Papyan et al.,, 2016; Papyan et al., 2017a, ) for neural networks in general, and CNN in particular.

To mitigate this problem for CNN-based imaging algorithms, we now present a physical interpretation of the CNN filters and feature maps in terms of the physics-based operators for seismic imaging. With such an understanding, we can use a physics-based rationale for the better design of CNN architectures in seismic migration and inversion.

Donoho, (2019) points out that “machine learning has a troubled relationship with understanding the foundation of its achievements well and its literature is admittedly corrupted by anti-intellectual and anti-scholarly tendencies”. Progress in advancing the capabilities of deep neural networks will be severely stymied unless its mathematical foundations are established. As a first step in this direction, Papyan et al., (2016) proposed that the forward modeling operation of CNN could be recast as finding the sparsest model ${\bf m}$ under the $L_{1}$ norm subject to honoring the data misfit constraint $||{{\bm{\Gamma}}{\bf m}-{\bf d}}||_{2}^{2}\leq\beta$:

[TABLE]

where ${\bf m}^{*}$ is the optimal solution of ${\bf m}$, ${\bm{\Gamma}}$ is the dictionary matrix, ${\bf d}$ represents the signal, and the scalar $\beta$ is the specified noise tolerance. The iterative solution to this problem is a series of forward-modeling operations of a neural network, where the mathematical operations of each layer consist two steps: a weighted summation of input values to give the vector ${\bf z}$ followed by a two-sided soft thresholding operation denoted as $\sigma({{\bf z}})$ (Papyan et al.,, 2016).

The sparse constraint problem defined in equation 1 is commonly seen in geophysics, for example, the least square migration (LSM) with a sparse constraint. LSM is an important seismic-imaging technique to produce images with better balanced amplitudes, fewer artifacts and better resolution than standard migration (Lailly,, 1983; Tarantola,, 1987; Schuster,, 1993; Nemeth et al.,, 1999; Chavent and Plessix,, 1999; Duquet et al.,, 2000; Feng and Schuster,, 2017; Schuster and Liu,, 2019). The sparse constraint is one of the important regularization terms used in solving the ill-conditioned least-squares problem (Sacchi and Ulrych,, 1995; De Roeck,, 2002; Kühl and Sacchi,, 2003; Wang and Sacchi,, 2005) and sparse LSM (SLSM) has been demonstrated to be effective for mitigation of aliasing artifacts and crosstalk in LSM (Wang and Sacchi,, 2007; Herrmann et al.,, 2009; Dutta,, 2017; Witte et al.,, 2017; Li and Gao,, 2018). The image-domain SLSM finds the sparsest reflectivity ${\bf m}^{*}$ with sparse constraints to minimize the objective function $||{{\bm{\Gamma}}{\bf m}-{\bf m}^{mig}}||_{2}^{2}$, where $\bm{\Gamma}$ is the Hessian matrix and ${\bf m}^{mig}$ is the migration image. Here, we can see SLSM shares the same problem as that defined in equation  1. Following the work of Papyan et al., (2016), we show that the sparse solution to the LSM problem reduces to the forward modeling operations of a multilayered neural network. The CNN filters and feature maps are shown to be analogous, but not equivalent, to the migration Green’s functions (Hessian) and the reflectivity distribution.

The standard SLSM algorithm needs to solve the wave equation, which is time-consuming. Motivated by the connection between sparse LSM and CNN, we propose the neural network version of sparse LSM, which does not need to solve the wave equation and is faster than standard SLSM. Instead of just finding the optimal reflectivity ${\bf m}^{*}$, we optimize for both the quasi-reflectivity $\bf m$ and the quasi-migration-Green’s functions $\bm{\Gamma}$. These quasi-migration-Green’s functions approximate the role of migration Green’s function (Schuster and Hu,, 2000) and are denoted as the convolutional filters in a convolutional neural network. As discussed in Appendix LABEL:A, the migration Green’s function is the point-scatterer response of the migration operator. The final image is denoted as the neural network least squares migration (NNLSM) estimate of the quasi-reflectivity distribution that honors the $L_{1}$ sparsity condition. The next section shows the connection between the multilayer neural network and the solution to the multilayer NNLSM problem. This is followed by the numerical examples with the synthetic models and field data from the North Sea.

2 Theory of Neural Network Least Squares Migration

The theory of standard image-domain LSM is first presented to establish the benchmark solution where the optimal reflectivity function minimizes the image misfit under the $L_{2}$ norm. This is then followed by the derivation of the sparse least squares migration (SLSM) solution for a single-layer network. The final two subsections derive the NNLSM solution for single-layer and multilayer networks, respectively.

2.1 Least Squares Migration

The least squares migration (LSM) problem can be defined (Schuster and Hu,, 2000; Schuster,, 2017) as finding the reflectivity coefficients $m_{i}$ in the $N\times 1$ vector ${\bf m}$ that minimize the $L_{2}$ objective function $\epsilon=1/2||{\bm{\Gamma}}{\bf m}-{\bf m}^{mig}||_{2}^{2}$,

[TABLE]

where ${\bm{\Gamma}}={\bf L}^{T}{\bf L}$ is the symmetric $N\times N$ Hessian matrix, ${\bf L}$ is the forward modeling operator, and ${\bf L}^{T}$ is the migration operator. Here, ${\bf m}^{mig}={{\bf L}^{T}}{\bf d}$ is the migration image computed by migrating the recorded data ${\bf d}$ with the migration operator ${\bf L}^{T}$. Alternatively, the image-domain LSM problem can also be defined as finding $\bf m$ that minimizes $\epsilon=1/2({\bf m}^{T}{\bm{\Gamma}}{\bf m}-{\bf m}^{T}{\bf m}^{mig})$, which has a more well-conditioned solution than the one in equation 2 (Schuster,, 2017). However, we will use equation 2 as the definition of the LSM problem in order to be consistent with the notation from Papyan et al., 2017a . The kernel associated with the Hessian matrix ${\bf L}^{T}{\bf L}$ is also known as the point scatterer response of the migration operator or the migration Green’s function (Schuster and Hu,, 2000). It is a square matrix that is assumed to be invertible, otherwise a regularization term is incorporated into the objective function.

A formal solution to equation 2 is

[TABLE]

where it is too expensive to directly compute the inverse Hessian ${\bm{\Gamma}}^{-1}$. Instead, a gradient method gives the iterative solution

[TABLE]

where $\alpha$ is the step length, $\bm{\Gamma}$ is symmetric, and ${\bf m}^{(k)}$ is the solution at the $k^{th}$ iteration. Typically, a regularization term is used to stabilize the solution, for example, the sparse constraint which will be introduced in the next subsection.

2.2 Sparse Least Squares Migration

The sparse least squares migration (SLSM) in the image domain is defined as finding the reflectivity coefficients $m_{i}$ in the $N\times 1$ vector ${\bf m}$ that minimize the objective function $\epsilon$ (Perez et al.,, 2013):

[TABLE]

where ${\bm{\Gamma}}={\bf L}^{T}{\bf L}$ represents the migration Green’s function (Schuster and Hu,, 2000), $\lambda>0$ is a positive scalar, ${\bf m}^{mig}={{\bf L}^{T}}{\bf d}$ is the migration image, and $S({\bf m})$ is a sparseness function. For example, the sparseness function might be $S({\bf m})=||{\bf m}||_{1}$ or $S({\bf m})=log(1+||{\bf m}||_{2}^{2})$.

The solution to equation 5 is

[TABLE]

which can be approximated by an iterative gradient descent method:

[TABLE]

Here, $S({\bf m})_{i}^{\prime}$ is the derivative of the sparseness function with respect to the model parameter $m_{i}$ and the step length is $\alpha$. Vectors and matrices are denoted by boldface lowercase and uppercase letters, respectively. When $S({\bf m})=||{\bf m}||_{1}$, the iterative solution in equation 7 can be recast as

[TABLE]

where, soft is the two-sided soft thresholding function (Papyan et al.,, 2016) derived in Appendix LABEL:B (see equation 22). Here, $\bm{\Gamma}={\bf L}^{T}{\bf L}$ is computed by solving the wave equation to get the forward modeled field and backpropagating the data by a numerical solution to the adjoint wave equation.

Equation 8 is similar to the forward modeling operation associated with the first layer of the neural network in Figure 1. That is, set $k=0$, ${\bf m}^{(0)}=0$, $\alpha=1$, and let the input vector be the scaled residual vector ${\bf r}=-({\bm{\Gamma}}{\bf m}^{(0)}-{\bf m}_{mig})={\bf m}_{mig}$ so that the first-iterate solution can be compactly represented by

[TABLE]

Here, the input vector ${\bf r}={\bf m}^{mig}$ is multiplied by the matrix ${\bm{\Gamma}}^{T}$ to give ${\bf z}={\bm{\Gamma}}^{T}{\bf r}$, and the elements of = ${\bf z}$ are then thresholded and shrunk to give the output ${{\bf m}}=\mbox{soft}({\bf z},\lambda)$. If we impose a positivity constraint for ${\bf z}$ and a shrinkage constraint so $\lambda$ is small, then the soft thresholding function becomes that of a one-sided threshold function, also known as the Rectified Linear Unit or ReLU function. To simplify the notation, the $\mbox{soft}({\bf z},\lambda)$ function or $\mbox{ReLu}({\bf z})$ function is replaced by $\sigma_{\lambda}({\bf z})$ so that equation 9 is given by

[TABLE]

For the ReLu function there is no shrinkage so $\lambda=0$. However, $\bm{\Gamma}$ in equation 10 is computed by forward and backward solutions to the wave equation. Unlike a neural network, the physics of wave propagation is included with the sparse LSM solution in equation 8.

2.3 Neural Network Least Squares Migration

We now propose the neural network version of SLSM that finds both $\bm{\Gamma}^{*}$ and ${\bf m}^{*}$ which minimize equation 5, which is equivalent to the convolutional sparse coding (CSC) problem. We denote the optimal solution $\bf m$ as the neural network least squares migration (NNLSM) image. Here, we assume that the migration image ${\bf m}^{mig}$ can be decomposed into components that have the form ${\bm{\Gamma}}_{1}{\bf m}_{1}$, where ${\bf m}_{1}$ represents a sparse quasi-reflectivity structure for the $1^{th}$ CNN layer in Figure 1 and ${\bm{\Gamma}}_{1}$ has a convolutional structure. The solution can be found by using the Alternating Direction Method of Multipliers (ADMM) method either in the Fourier domain (Heide et al.,, 2015) or in the space domain (Papyan et al., 2017b, ), which alternates between finding $\bm{\Gamma}^{*}$ (dictionary learning problem) and then finding ${\bf m}^{*}$ (sparse pursuit problem). The key difference between least squares sparse inversion and NNLSM is that $\bm{\Gamma}$ is not computed by numerical solutions to the wave equation. Instead, the coefficients of $\bm{\Gamma}$ are computed by the usual gradient-descent learning algorithm of a convolutional neural network. For this reason, we denote $\bf m^{*}$ as the quasi-reflectivity distribution and $\bm{\Gamma}^{*}$ as the quasi-migration Green’s function.

Appendix LABEL:C shows the general solution for NNLSM for a single-layer neural network, where the optimal $\bm{\Gamma}^{*}$ is composed of the quasi-migration Green’s functions, which are denoted as convolutional filters in the machine learning terminology (Liu and Schuster,, 2018, 2019). Each filter is used to compute a feature map that corresponds to a sub-image of quasi-reflection coefficients in the context of LSM.

We now compute the NNLSM image for a 1D model, where we assume ${\bf m}^{mig}$ is a N-dimensional vector which can be expressed as,

[TABLE]

Here, ${\bm{\gamma}}_{i}$ is the $i^{th}$ local filter with length of $n_{0}$, ${\bf m}^{\prime}_{1i}$ is the $i^{th}$ feature map, “*” denotes the convolution operator and $k_{0}$ is the number of the filters. Alternatively, following Figure 2a, equation 11 can be written in matrix form as ${\bf m}^{mig}={\bm{\Gamma}}_{1}{\bf m}_{1}={\bm{\Gamma}}_{1}^{\prime}{\bf m}_{1}^{\prime}$ (Papyan et al., 2017a, ), where ${\bm{\Gamma}}_{1}$ is a convolutional matrix containing in its columns the $k_{0}$ filters with all of their shifts. ${\bm{\Gamma}}_{1}^{\prime}$ is a concatenation of banded and circulant 222We shall assume throughout this paper that boundaries are treated by a periodic continuation, which gives rise to the cyclic structure. matrices, which is the same as ${\bm{\Gamma}}_{1}$ except that the order of the columns is different. ${\bf m}_{1}^{\prime}$ is a concatenation of the feature map vectors ${\bf m}_{1i}^{\prime}$ for $i=1,2,\cdots,k_{0}$.

The advantage of NNLSM is that only inexpensive matrix-vector multiplications are used and no expensive solutions to the wave equation are needed for backward and forward wavefield propagation. As will be seen later, convolutional filters that appear to be coherent noise can be excluded for denoising the migration image.

2.4 Multilayer Neural Network LSM

The multilayer NNLSM is a natural extension of the single-layer NNLSM. For NNLSM, the migration image ${\bf m}^{mig}$ can be expressed as ${\bf m}^{mig}={\bm{\Gamma}}_{1}{\bf m}_{1}$ (Figure 2a), where there are $k_{0}$ filters in ${\bm{\Gamma}}_{1}$ and $k_{0}$ sub-quasi-reflectivity images in ${\bf m}_{1}$. Following Sulam et al., (2018), we can cascade this model by imposing a similar assumption to the sparse representation ${\bf m}_{1}$, i.e., ${\bf m}_{1}={\bm{\Gamma}}_{2}{\bf m}_{2}$, for a corresponding convolutional matrix ${\bm{\Gamma}}_{2}$ with $k_{1}$ local filters and a sparse sub-quasi-reflectivity image ${\bf m}_{2}$, as depicted in Figure 2b. In this case, the filter size is $n_{1}\times k_{0}$ and there are $k_{1}$ sub-quasi-reflectivity images in ${\bf m}_{2}$.

Similar to the derivation by Papyan et al., 2017a and Sulam et al., (2018), the multilayer neural network LSM problem is defined as the following.

[TABLE]

where ${\bm{\Gamma}}_{i}$ is the $ith$ Hessian matrix in the $ith$ layer. The first iterate solution to the above system of equations can be cast in a form similar to equation 10, except we have

[TABLE]

which is a repeated concatenation of the two operations of a multilayered neural network: matrix-vector multiplication followed by a thresholding operation. In all cases, we use a convolutional neural network where different filters are applied to the input from the previous layer to give feature maps associated with the next layer, as shown in Figure 1.

For a perfect prediction of the migration image, ${\bf m}^{mig}$ can also be approximated as ${\bf m}^{mig}={\bm{\Gamma}}_{1}{\bm{\Gamma}}_{2}\dots{\bm{\Gamma}}_{N}{\bf m}_{N}$. We refer to ${\bm{\Gamma}}^{(i)}$ as the effective filter at the $i^{th}$ level,

[TABLE]

so that

[TABLE]

The next section tests the effectiveness of NNLSM on both synthetic data and field data.

3 Numerical Results

We now present numerical simulations of NNLSM. Instead of only determining the optimal reflectivity $\bf m$ as computed by SLSM, the NNLSM method computes both quasi-reflectivity $\bf m$ and the elements of the Hessian matrix ${\bm{\Gamma}}=\bf L^{T}{\bf L}$. Each block of ${\bm{\Gamma}}$ is considered to be the segment response function (SSF) of the migration operator rather than the point spread function (PSF). If the actual Green’s functions are used to construct ${\bm{\Gamma}}$ then each column of the Hessian matrix is the point scatterer response of the migration operator (Schuster and Hu,, 2000). In contrast, the NNLSM Hessian is composed of blocks, where each block is the segment scatterer response of the migration operator. An example will be shown later where a segment of the reflector is migrated to give the migration segment response of the migration operator. The computational cost for computing SSF’s is several orders of magnitude less than that for PSFs because no solutions to the wave equation are needed. The penalty, however, is that the resulting solution $\bf m$ is not the true reflectivity, but a sparse representation of it we denote as the quasi-reflectivity distribution.

Using the terminology of neural networks, we can also denote the sparse sub-quasi-reflectivity images as feature maps. Each block in $\bm{\Gamma}$ will be denoted as a filter. Therefore the vector output of ${\bm{\Gamma}}\bf m$ can be interpreted as a sum of filter vectors ${\bm{\gamma}}_{i}$ weighted by the coefficients in $\bf m$, where ${\bm{\gamma}}_{i}$ is the $i^{th}$ column vector of ${\bm{\Gamma}}$.

3.1 Three-layer Velocity Model

The interpretation of feature maps and filters can be understood by computing them for the Figure 3a model. The grid size of the model is 101$\times$101, and the grid interval is 10 m in both the x and z directions. There are 26 shots evenly spaced at a distance of 40 m on the surface, and each shot is recorded by 101 receivers with a sampling interval of 10 m. Figure 3b show the reverse time migration (RTM) image.

The first test is for a 1D model where we extract the image located at X$=0.5$ km, which is displayed as the blue curve in Figure 3c. The red curve in Figure 3c is the reflectivity model. Assume that there is only one filter in ${\bm{\Gamma}}$ and it extends over the depth of 400 m (41 grid points). We now compute the NNLSM image by finding the optimal ${\bf m}$ and ${\bm{\Gamma}}$ by the two-step iterative procedure denoted as the alternating descent method (see Liu and Schuster, (2018) and Liu et al., (2018)). The computed filter $\bm{\gamma}$i is shown in Figure 4a where the phase of the filter ${\bm{\gamma}}_{i}$ is nonzero. If we use a filter with a non-zero time lag to calculate its feature map ${\bf m}$, the phases of the feature map and the true reflectivity ${\bf m}$ will be different. So, we need to modify the time lag and polarity of the basis function $\tilde{\bm{\Gamma}}_{i}$. The modified basis function is shown in Figure 4b, and its coefficients are displayed as the blue curve in Figure 4c. Compared with the true reflectivity ${\bf m}$ (red curve in Figure 4), the feature map can give the correct positions but also give the wrong values of the reflectivity distribution.

Next, we perform a 2D test where the input is the 2D migration image in Figure 3b. Three 35-by-35 (grid point) filters are learned (see Figure 5a). The modified filters are shown in Figure 5b. Appendix LABEL:D describes how we align the filters by using the cross-correlation method. The feature maps of these three filters are displayed in Figures 6a-6c. Figure 6d shows the sum of these three feature maps. It is evident that the stacked feature maps can estimate the correct locations of the reflectivity spikes.

3.2 SEG/EAGE Salt Model

The multilayer NNLSM procedure (see equation 12) is now applied to the migration image associated with the 2D SEG/EAGE salt velocity model in Figure 7a. The grid size of the model is 101 grid points in both the z- and x-directions. The grid interval is 40 m in the x-direction and 20 m in the z-direction. Figure 7b shows the reverse time migration (RTM) image. The multilayer NNLSM consists of three convolutional layers: the first one contains 15 basis functions, i.e., filters, of size 11$\times$11 grid points, the second one consists of 15 basis functions with dimensions 11$\times$11$\times$15, and the last one contains contains 15 basis function of dimensions 11$\times$11$\times$15. Equation 12 is solved for both ${{\bf m}_{i}}$ and ${\bm{\Gamma}}_{i}$ ($i\in 1,2,3$) by the two-step iterative procedure denoted as the alternating descent method. The multilayered structure is shown in Figure 8, where the black dots in ${\bf m}_{i}$ represent the nonzero values of the quasi-reflectivity distribution. The effective basis functions computed for these layers are shown in Figures 7c-7e, where the yellow, red and green boxes indicate the sizes of the effective basis functions, which can be considered as quasi-migration Green’s functions. It indicates that the basis functions of the first layer ${\bm{\Gamma}}_{1}$ contain very simple small-dimensional edges, which are called “atoms” by Sulam et al., (2018). The non-zeros of the second group of basis functions ${\bm{\Gamma}}_{2}$ combine a few atoms from ${\bm{\Gamma}}_{1}$ to create slightly more complex edges, junctions and corners in the effective basis functions in ${\bm{\Gamma}}^{(2)}$. Lastly, ${\bm{\Gamma}}_{3}$ combines atoms from ${\bm{\Gamma}}^{(2)}$ in order to reconstruct the more complex parts of the migration image. The corresponding stacked coefficient images, also known as feature maps, are shown in Figures 7f-7h, which give the quasi-reflectivity distributions. The reconstructed migration images are shown in Figures 7i-7k.

For comparison, we computed the standard LSM image using the deblurring method described in Chen et al., (2017, 2019). Here, the deblurring filter size is 17x17 grid points (black boxes in Figure 9) and computed for a 50x50 grid (red boxes in Figure 9) of evenly spaced point scatterers with the same migration velocity model as used for the data migration in Figure 7a. The standard LSM images for the first and 50th iterations are shown in Figures  10b and  10c, respectively, next to the NNLSM image in Figure 10d. It is clear that the NNLSM image is better resolved than the LSM image, although there are discontinuities in some of the NNLSM interfaces not seen in the LSM image. Some of the detailed geology is lost in the LSM image as seen in the wiggly interface in the red-rectangular area of Figure 10. The practical application of the NNLSM image is that it might serve a super-resolved attribute image that can be combined with other attributes to delineate geology. For example, combining the depth-slice of the NNLSM image with a spectral decomposition image (Aarre,, 2016) can help delineate the lithological edges of meandering channels.

NNLSM can filter out random and coherent noises in the migration image after reconstructing the migration image by eliminating the noisy learned basis functions and their coefficients in the NNLSM image. For example, Figure 11a shows the RTM image with a sparse acquisition geometry so that the image contains a strong acquisition footprint. The reconstructed migration image in Figure 11b shows significant mitigation of this noise in Figure 11a. However, the migration swing noise is still prominent near the red arrows in Figure 11b. Such noise is reconstructed from the noisy basis function shown in Figure 12a and the coefficients in Figure 12b. Figure 12c is the image reconstructed by correlating the basis function in Figure 12a with the coefficients in Figure 12b. After filtering out the basis functions from noise, the reconstructed image is shown in Figure 11c, which is free from aliasing noise at the locations indicated by the red arrows.

3.3 North Sea Data

We apply the NNLSM method to field data collected in the North Sea (Schroot and Scüttenhelm,, 2003), where the time migration image is shown in Figure 13a. The time axis is gridded with 213 evenly-spaced points and there are 301 grid points along the x-axis. We compute 21 13-by-5 (grid point) convolutional basis functions, i.e. filters $\bm{\gamma}_{i}$ for ($i=1,2,...21$), by the NNLSM procedure (see Figure 13b). These filters approximate the dip-filtered migration Green’s functions, and the basis function is marked as the yellow boxes in Figure 13a and 13b. The stacked feature maps (quasi-reflectivity distribution) are displayed in Figure 13c. It is evident that the stacked feature maps can provide a high-resolution migration image. After reconstruction from the learned filters and feature maps, the migration image is shown in Figure 13d with less noise.

Finally, we apply NNLSM to a time slice of the migration image, which is shown in Figure 14a, and the image size is 301 by 301 gridpoints. Figure 14b shows the 21 13$\times$5 filters estimated by the NNLSM procedure. The stacked feature map is displayed in Figure 14c, which may be used as a superresolution attribute image for high-resolution delineation of geologic bodies. The reconstructed migration image is shown in Figure 14d and we can see there is less noise.

4 Discussion

The forward modeling for a multilayered neural network is shown to be equivalent to a single-iterate solution of a multilayered LSM problem. This assumes positivity and shrinkage constraints on the soft thresholding operation, so it reduces to the ReLu operation. This equivalence relates the physics of seismic imaging to architectural features in the neural network.

• •

The size of the filters in the first layer should be about the same size as the Green’s function for that model. Experiments with numerical models suggest that this size is approximately one-to-two wavelengths. In this case, the filter is interpreted as an approximation to the migration Green’s function, except it is that for a reflecting segment. Thus, we interpret the approximate migration Green’s function as a migration segment spread function (SSF) rather than a migration point spread function. Sulam et al., (2018) classifies each feature in the first layer as an atom which takes on the role of a SSF.

• •

The output of the first layer provides the small scale, i.e., high-wavenumber, features associated with the input data. For an input migration image, the feature maps of the first layer resemble sub-quasi-reflectivity maps of the subsurface. Adding the sub-quasi-reflectivity maps together gives a close approximation to the actual reflectivity model as shown in Figures 6d and 7f.

• •

The output of the second layer is a weighted sum of the first-layer features, which create sparser feature maps. Sulam et al., (2018) classifies the concatenation of the filters from the first and second layers as molecules (see equation 14). In the migration problem, the resulting filters are SSFs for even larger segments of the original reflector boundaries. The feature maps of the third layer are a weighted sum of the second layer’s features to produce even the sparsest feature maps. For migration, the final feature maps are very sparse while the concatenated filters are associated with large-scale features of the migration image.

• •

The computational cost of computing NNLSM images is significantly less than that for LSM images because no solutions of the wave equation are needed. For example, we consider the 2D FDTD forward modeling of the acoustic equation with an eighth-order scheme in space and a second-order scheme in time, and its computational complexity is $O(n_{t}n^{2})$ for one shot, where $n$ is the number of grid points in one direction and $n_{t}$ is the number of the time steps. According to Plessix, (2007), $n_{t}$ is approximate $30\ n$ to satisfy the stability condition and to make sure there is enough recoding time when $v_{max}/v_{min}\simeq=3$, where $v_{max}$ and $v_{min}$ are the maximum and minimum velocities, respectively. So, The complexity of 2D FDTD forward modeling of the acoustic equation is $O(n_{s}n^{3})$ where $n_{s}$ is the number of the shots. The complexity of LSRTM is $O(6n_{s}N^{iter}n^{3})$ (Schuster,, 2017), where $N^{iter}$ is the iteration number. For NNLSM, the complexity is $O(N^{iter}n^{2}\log n)$ according to Heide et al., (2015) and can be reduced to $O(N^{iter}n^{2})$ if a local block coordinate-descent algorithm is used (Zisselman et al.,, 2019).

The 1D NNLSM can be interpreted as a blind deconvolution (BD) problem in seismic data processing (Kaaresen and Taxt,, 1998; Bharadwaj et al.,, 2018). It can be seen in Figure 5 that the filter of NNLSM is the source wavelet of BD and the coefficients of NNLSM are the quasi-reflectivity coefficients of BD. However, NNLSM can have more than one filter and the filters can be 2D or 3D filters. We show the relationship between BD, sparse LSM, NNLSM, and CNN in Figure 15.

NNLSM is an unsupervised learning method. Compared to a supervised learning method, it does not heavily depend on the huge amount of training data. However, it may need human intervention for inspecting the quality of the quasi-reflectivity images and we need to align the filters to get more consistent quasi-reflectivity.

In Figure 14, we apply 2D NNLSM to a time slice of the 3D North Sea Data. The 2D NNLSM method can be extended to a 3D implementation to learn a set of 3D filters. Incorporating the third dimension of information from 3D filters will likely lead to a better continuity of the reflectors in Figure 14c. However, the computational cost will increase by more than a multiplicative factor of $n$.

5 Conclusion

Neural network least squares migration finds the optimal quasi-reflectivity distribution and quasi-migration-Green’s functions that minimize a sum of migration misfit and sparsity regularization functions. The advantages of NNLSM over standard LSM are that its computational cost is significantly less than that for LSM and it can be used for filtering both coherent and incoherent noise in migration images. A practical application of the NNLSM image is as an attribute map that provides superresolution imaging of the layer interfaces. This attribute image can be combined with other attributes to delineate both structure and lithology in depth/time slices of migration images. Its disadvantage is that the NNLSM quasi-reflectivity image is only an approximation to the actual reflectivity distribution.

A significant contribution of our work is that we show that the filters and feature maps of a multilayered CNN are analogous to the migration Green’s functions and reflectivity distributions. For the first time we now have a physical interpretation of the filters and feature maps in deep CNN in terms of the operators for seismic imaging. Such an understanding has the potential to lead to better architecture design of the network and extend its application to waveform inversion. In answer to Donoho’s plea for more rigor, NNLSM represents a step forward in establishing the mathematical foundation of CNN in the context of least squares migration.

6 Acknowledgements

The research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. We are grateful to the sponsors of the Center for Subsurface Imaging and Modeling Consortium for their financial support. For computer time, this research used the resources of the Supercomputing Laboratory at KAUST and the IT Research Computing Group. We thank them for providing the computational resources required for carrying out this work.

\append

Migration Green’s Function

Schuster and Hu, (2000) show that the poststack migration (Yilmaz,, 2001) image $m({\bf x})^{mig}$ in the frequency domain is computed by weighting each reflectivity value $m({\bf z})$ by $\Gamma({\bf x}|{\bf z})$ and integrating over the model-space coordinates ${\bf z}$:

[TABLE]

where $\eta$ represents terms such as the frequency variable raised to the 4th power. The migration Green’s function $\Gamma({\bf x}|{\bf z})$ is given by

[TABLE]

Here we implicitly assume a normalized source wavelet in the frequency domain, and $D_{model}$ and $D_{data}$ represent the sets of coordinates in, respectively, the model and data spaces. The term $G({\bf x}^{\prime}|{\bf x})=e^{i\omega\tau_{xx^{\prime}}}/||{\bf x}-{\bf x}^{\prime}||$ is the Green’s function for a source at ${\bf x}$ and a receiver at ${\bf x}^{\prime}$ in a smoothly varying medium333If the source and receiver are coincident at ${\bf x}$ then the zero-offset trace is represented by the squared Green’s function $G({\bf x}|{\bf x}^{\prime})^{2}$.. The traveltime $\tau_{xx^{\prime}}$ is for a direct arrival to propagate from ${\bf x}$ to ${\bf x}^{\prime}$.

The physical interpretation of the kernel $\Gamma({\bf x}^{\prime}|{\bf x})$ is that it is the migration operator’s444This assumes that the zero-offset trace is generated with an impulsive point source with a smoothly varying background velocity model, and then migrated by a poststack migration operation. It is always assumed that the direct arrival is muted and there are no multiples. response at ${\bf x}^{\prime}$ to a point scatterer at ${\bf x}$, otherwise known as the MGF or the migration Green’s function (Schuster and Hu,, 2000). It is analogous to the point spread function (PSF) of an optical lens for a point light source at ${\bf x}$ in front of the lens and its optical image at ${\bf x}^{\prime}$ behind the lens on the image plane. In discrete form, the modeling term $[{\bm{\Gamma}{\bf m}}]_{i}$ in equation 16 can be expressed as

[TABLE]

with the physical interpretation that $[{\bm{\Gamma}}{\bf m}]_{i}$ is the migration Green’s function response at ${\bf x}_{i}$. An alternative interpretation is that $[{\bm{\Gamma}}{\bf m}]_{i}$ is the weighted sum of basis functions $\Gamma({\bf x}_{i}|{\bf z}_{j})$ where the weights are the reflection coefficients $m_{j}$ and the summation is over the $j$ index. We will now consider this last interpretation and redefine the problem as finding both the weights $m_{i}$ and the basis functions $\Gamma({\bf x}_{i}|{\bf z}_{j})$. This will be shown to be equivalent to the problem of a fully connected (FCN) neural network.

\append

Soft Thresholding Function Define the sparse inversion problem as finding the optimal value $x^{*}$ that minimizes the objective function

[TABLE]

where the $L_{1}$ norm demands sparsity in the solution ${\bf x}$. An example is where ${\bf z}$ is a noisy $M\times N$ image such that ${\bf z}={\bf x}+\mbox{noise}$, and we seek the optimal vector ${\bf x}$ that satisfies equation 19. Here, the noisy $M\times N$ image has been flattened into the tall $MN\times 1$ vector ${\bf z}$.

The stationary condition for equation 19 is

[TABLE]

where

[TABLE]

Equations 20-21 can be combined to give the optimal $x^{*}$ expressed as the two-sided ReLu function

[TABLE]

More generally, the iterative-soft-threshold-algorithm (ISTA) that finds ${\bf x}^{*}$

[TABLE]

is

[TABLE]

There are several more recently developed algorithms that have faster convergence properties than ISTA. For example, FISTA (Fast-ISTA) has quadratic convergence (Beck and Teboulle,, 2009). \appendNeural Network Least Squares Migration

The neural network least squares migration (NNLSM) algorithm in the image domain is defined as solving for both the basis functions $\tilde{\Gamma}({\bf x}_{i}|{\bf x}_{j})$ and $\tilde{m}_{j}$ that minimize the objective function defined in equation 5. In contrast, SLSM only finds the least squares migration image in the image domain and uses the pre-computed migration Green’s functions that solve the wave equation.

The NNLSM solution is defined as

[TABLE]

where now both $\tilde{\bm{\Gamma}}^{*}$ and $\tilde{\bf m}^{*}$ are to be found. The functions with tildes are mathematical constructs that are not necessarily identical to those based on the physics of wave propagation as in equation 5.

The explicit matrix-vector form of the objective function in equation 25 is given by

[TABLE]

and its Fréchet derivative with respect to $\tilde{\Gamma}({\bf x}_{i^{\prime}}|{\bf z}_{j^{\prime}})$ is given by

[TABLE]

The iterative solution of equation 25 is given in two steps (Olshausen and Field,, 1996).

Iteratively estimate $\tilde{m}_{i}$ by the gradient descent formula used with SLSM:

[TABLE]

However, one migration image ${\bf m}^{mig}$ is insufficient to find so many unknowns. In this case the original migration image is broken up into many small pieces so that there are many migration images to form examples from a large training set. For prestack migration, there will be many examples of prestack migration images, one for each shot, and the compressive sensing technique denoted as VISTA (Ahmad et al.,, 2015) is used for the calculations. 2. 2.

Update the basis functions $\tilde{\Gamma}({\bf x}_{i}|{\bf z}_{j})$ by inserting equation 27 into the gradient descent formula to get

[TABLE]

It is tempting to think of $\tilde{\Gamma}({\bf x}|{\bf x}^{\prime})$ as the migration Green’s function and $\tilde{m}_{i}$ as the component of reflectivity. However, there is yet no justification to submit to this temptation and so we must consider, unlike in the SLSM algorithm, that $\tilde{\Gamma}({\bf x}|{\bf x}^{\prime})$ is a sparse basis function and $\tilde{m}_{i}$ is its coefficient. To get the true reflectivity then we should equate equation 18 to $\sum_{j}\tilde{\Gamma}({\bf x}_{i},{\bf x}_{j})\tilde{m}_{j}$ and solve for $m_{j}$.

\append

Alignment of the Filters

To align the learned filters, we first choose a “target” filter, which is denoted as a 2D matrix ${\bf A}$ with the size of $M\times N$. Then we try to align all the other filters with the target filter through their cross-correlation. For example, if we choose one filter denoted as matrix ${\bf B}$ with the same size as ${\bf A}$, we can get the cross-correlation matrix $\bf C$ with its elements defined as,

[TABLE]

where $-M<i<M$ and $-N<j<N$. $a_{i,j}$, $b_{i,j}$ and $C_{i,j}$ indicate the element at position $(i,j)$ in matrices $\bf A$, $\bf B$ and $\bf C$, respectively. The location of the maximum absolute value of the elements in matrix $\bf C$ indicates how much should we shift filter $\bf B$ to filter $\bf A$ in each direction. Figure 16 shows the calculation of the cross-correlation matrix $\bf C$ for two filters $A$ and $B$. $c_{1,0}$ (or $C_{4,3}$) is the maximum absolute value of the elements in matrix $\bf C$, which indicates filter $\bf B$ should be shifted 1 position along the first direction. Here, we need to pad zeros along all the dimensions of filter $\bf B$ before shifting it, which is displayed in Figure 17.

Figure 18a shows the learned filters with a size of $17\times 9$ from the migration image of the SEG/EAGE salt model. Filter No. 7 (yellow box) is chosen as the target filter. The aligned filters are shown in Figure 18b without zero padding and the stacked feature maps from the original and aligned filters are displayed in Figures 18c and 18d, respectively. It is evident that the reflector interfaces from the aligned filters are more continuous especially in the red box compared with those of the original filters.