DuRIN: A Deep-unfolded Sparse Seismic Reflectivity Inversion Network

We broadly classify the prior work related to reflectivity inversion into two categories, namely, optimization-based techniques and data-driven techniques.

The limitations of $\ell_{1}$ regularization [20, 21, 22], in addition to the challenge of estimating the sparsity of seismic reflection data [8, 19], can be overcome through nonconvex regularization [19]. However, the corresponding nonconvex cost suffers from local minima, and optimization becomes quite challenging. The hyper-parameters of most iterative techniques are set heuristically, which is likely to yield suboptimal solutions. Machine learning approaches outperform conventional techniques in reflectivity inversion [3], especially in support recovery [34], which is given higher priority [5]. Deep-unrolled architectures [51], which belong to a class of model-based neural networks, are not entirely dissociated from the underlying physics of the problem, which is a potential risk in elementary data-driven approaches such as feedforward neural networks [3, 34, 27]. These factors motivate us to employ a data-driven nonconvex regularization strategy and solve the reflectivity inversion problem through deep-unrolled architectures.

The contributions of this paper are summarized as follows.

1. 1.

   We construct an optimization cost for sparse seismic reflectivity inversion based on the weighted counterpart of the minimax-concave penalty (MCP) [24], where each component of the sparse reflectivity ${\boldsymbol{x}}$ is associated with a different weight.

2. 2.

   The resulting optimization algorithm, the iterative firm-thresholding algorithm (IFTA) [58] is unfolded into the deep-unfolded reflectivity inversion network (DuRIN). To the best of our knowledge, model-based architectures have not been explored for solving the seismic reflectivity inversion problem. In fact, deep-unrolling has not been employed for solving seismic inverse problems [28].

3. 3.

   The efficacy of the proposed formulation with respect to the state of the art is demonstrated over synthetic $1$-D and $2$-D data and on simulated as well as real datasets.

Formulating the optimization problem using the weighted MCP regularizer defined above leads to the objective function:

The objective function ${J}$, data-fidelity term $F$, and the regularizer $g$ satisfy the following properties, which are essential for minimizing ${J}$:

1. P1.

   $F:\mathbb{R}^{n}\rightarrow\left(-\infty,\infty\right]$ is proper, closed, and $L$-smooth, i.e., $\|\nabla F({\boldsymbol{x}})-\nabla F({\boldsymbol{y}})\|_{2}\leq L\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{2}$.

2. P2.

   $g$ is lower semi-continuous.

3. P3.

   ${J}$ is bounded from below, i.e., inf ${J}>-\infty$.

Next, we optimize the problem stated in (7) through the Majorization-Minimization (MM) approach [66], and the resulting algorithm, the iterative firm-thresholding algorithm (IFTA) [58]. Unfolding the iterations of IFTA results in a learnable network called deep-unfolded reflectivity inversion network (DuRIN), which has a similar architecture to FirmNet [58].

Due to P1, there exists $\eta<1/L$ such that $F({\boldsymbol{x}})$ is upper-bounded locally by a quadratic expansion about ${\boldsymbol{x}}={\boldsymbol{x}}^{(k)}$ as:

where ${Q^{(k)}\left(\boldsymbol{x}\right)}=F({\boldsymbol{x}}^{(k)})+({\boldsymbol{x}}-{\boldsymbol{x}}^{(k)})^{\intercal}\nabla F({\boldsymbol{x}}^{(k)})+\frac{1}{2\eta}{\left\|{\boldsymbol{x}}-{\boldsymbol{x}}^{(k)}\right\|}_{2}^{2}$. The majorizer to the objective function $J$ at ${\boldsymbol{x}}^{(k)}$ is given by:

such that $H^{(k)}({\boldsymbol{x}})\geq J(\boldsymbol{x})$. The update equation for minimizing $H^{(k)}({\boldsymbol{x}})$ to obtain ${\boldsymbol{x}}^{(k+1)}$ is given by,

The proximal operator of the penalty function $g_{{\gamma}_{i}}({x}_{i};{\mu}_{i})$ is the firm-thresholding operator [67] given by

The update step at the $(k+1)^{th}$ iteration is given by

where ${\boldsymbol{h}}^{\prime}$ is the flipped version of ${\boldsymbol{h}}$. Algorithm (1) lists the steps involved in solving the optimization problem (7). The update in (12) involves convolutions followed by nonlinear activation (the firm threshold in this case), and can therefore be represented as a layer in a neural network.

As mentioned in the previous section, the update in IFTA (12) can be interpreted as a layer in a neural network, and hence, we unroll the iterations of the algorithm into layers of a neural network to solve the problem given in (7). The proposed deep-unrolled architecture for sparse seismic reflectivity inversion is called Deep-unfolded Reflectivity Inversion Network (DuRIN), which resembles FirmNet [58]. The input-output for each layer in DuRIN is given by

where $\mathbf{W}=\eta\mathop{{\boldsymbol{h}}^{\prime}}$ and $\mathbf{S}=\mathbf{I}-\eta\mathop{{\boldsymbol{h}}^{\prime}{\boldsymbol{h}}}$ [50] are initialized as Toeplitz matrices, and are dense and unstructured in the learning stage. The parameters that need to be learned are the matrices $\mathbf{W}$ and $\mathbf{S}$ and the thresholds (${\boldsymbol{\mu}},{\boldsymbol{\gamma}}$), subject to ${\boldsymbol{\mu}}>{\boldsymbol{0}},{\boldsymbol{\gamma}}>{\boldsymbol{1}}$, where ${\boldsymbol{0}}$ is the null vector and ${\boldsymbol{1}}$ is the vector of all ones. For training, we employ the smooth $\ell_{1}$ loss $\left(\text{Smooth}~{}\ell_{1}=\beta\left\|{{\boldsymbol{x}}-\hat{{\boldsymbol{x}}}}\right\|_{1}+(1-\beta)\left\|{{\boldsymbol{x}}-\hat{{\boldsymbol{x}}}}\right\|_{2}^{2},~{}0<\beta\leq 1\right)$ as the training cost, computed between the true reflectivity ${\boldsymbol{x}}$ and prediction $\hat{{\boldsymbol{x}}}$. In our experimental setup, $\beta=1$, and it is not a trainable parameter.

To improve amplitude recovery of DuRIN during inference, we re-estimate the amplitudes of the estimated sparse vector $\hat{{\boldsymbol{x}}}$ over the supports given by DuRIN as: $\hat{{\boldsymbol{x}}}_{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}}=\mathbf{H}_{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}}^{\dagger}{\boldsymbol{y}}$, where $\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}$ is the support, or the non-zero locations, of $\hat{\boldsymbol{x}}$, and $\mathbf{H}_{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}}^{\dagger}$ is the pseudo-inverse of the Toeplitz matrix $\mathbf{H}$ of the kernel ${\boldsymbol{h}}$ over the support $\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}$. The steps of DuRIN for $k_{max}$ layers are listed in Algorithm 2. Figure 4 illustrates a 3-layer architecture for DuRIN.

We consider two variants of DuRIN, depending on the value of $\gamma_{i},~{}\forall~{}i\in[\![1,n]\!]$. The generalized variant of DuRIN based on IFTA [58], where ${\gamma}_{i}$ are different across components and layers, is referred to as DuRIN-$1$, and corresponds to FirmNet [58]. As $\gamma_{i}\to\infty$, the firm-thresholding function approaches the soft thresholding function [22]. We call this variant DuRIN-$2$, which resembles LISTA [50].

We evaluate the results using statistical parameters that measure the amplitude and support recovery between the ground-truth sparse vector ${\boldsymbol{x}}$ and the predicted sparse vector $\hat{{\boldsymbol{x}}}$.

1. 1.

   Correlation Coefficient (CC) [71]:

   $\displaystyle\text{CC}=\frac{\langle{\boldsymbol{x}},\hat{{\boldsymbol{x}}}\rangle-\langle{\boldsymbol{x}},{{\boldsymbol{1}}}\rangle\langle\hat{{\boldsymbol{x}}},{{\boldsymbol{1}}}\rangle}{\sqrt{\left(\left\|{{\boldsymbol{x}}}\right\|_{2}^{2}-{\langle{\boldsymbol{x}},{{\boldsymbol{1}}}\rangle}^{2}\right)\left(\left\|\hat{{\boldsymbol{x}}}\right\|_{2}^{2}-{\langle\hat{{\boldsymbol{x}}},{{\boldsymbol{1}}}\rangle}^{2}\right)}},$

   where $\langle\cdot,\cdot\rangle$ denotes inner product, and ${\boldsymbol{1}}$ denotes a vector of all ones.

2. 2.

   Relative Reconstruction Error (RRE) and Signal-to-Reconstruction Error Ratio (SRER), defined as follows:

   $\displaystyle\text{RRE}=\frac{{\left\|\hat{{\boldsymbol{x}}}-{\boldsymbol{x}}\right\|}_{2}^{2}}{{\left\|{\boldsymbol{x}}\right\|}_{2}^{2}},\quad\text{SRER}=10\mathop{\log\limits_{10}\left(\frac{{\left\|{{\boldsymbol{x}}}\right\|_{2}^{2}}}{{\left\|{\hat{{\boldsymbol{x}}}-{{\boldsymbol{x}}}}\right\|_{2}^{2}}}\right){\text{dB}}}.$

3. 3.

   Probability of Error in Support (PES):

   $\displaystyle\text{PES}=\left(\frac{{{\text{max}}(\left|{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}}\right|,\left|{{\mathcal{S}}_{\boldsymbol{x}}}\right|)-\left|{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}\mathop{\cap}{\mathcal{S}_{\boldsymbol{x}}}}\right|}}{{{\text{max}}(\left|{\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}}\right|,\left|{\mathcal{S}_{\boldsymbol{x}}}\right|)}}\right),$

   where $|\cdot|$ denotes the cardinality of the argument, and ${\mathcal{S}}_{\boldsymbol{x}}$ and $\hat{\mathcal{S}}_{\hat{\boldsymbol{x}}}$ denote the supports of ${\boldsymbol{x}}$ and $\hat{{\boldsymbol{x}}}$, respectively.

The generation or acquisition of training data appropriately representing observed seismic data is crucial for the application of a data-driven approach to the reflectivity inversion problem [34]. We generate synthetic training data as $5\times 10^{5}$ seismic traces, each of length $300$ samples, obtained by the convolution between $1$-D reflectivity profiles and a Ricker wavelet. The reflectivity profiles consist of $200$ samples each, and are padded with zeros before convolution such that their length is equal to that of the seismic traces [3]. Amplitudes of the reflection coefficients range from $-1.0$ to $1.0$, and the sparsity factor ($k$), i.e., the ratio of the number of non-zero elements to the total number of elements in a trace, is set to $0.05$. The locations of the reflection coefficients are picked uniformly at random, with no constraint on the minimum spacing between two reflectors/spikes, or in other words, the reflectors can be away only by one sampling interval. The reflectivity values, chosen randomly from the defined range of $-1.0$ to $1.0$, are assigned to the selected locations. The sampling interval, the increment in the amplitudes of the reflection coefficients, and the wavelet frequency vary based on the dataset.

The initial hyperparameters and the optimum number of layers in the network also vary depending on the parameters of the dataset, such as the sampling interval and the wavelet frequency. We initially consider six layers for both DuRIN-$1$ and DuRIN-$2$ for training and keep appending five layers at a time to arrive at the optimum number of layers. For the training phase, we set the batch size to $200$, use the ADAM optimizer [72] and set the learning rate to $1\times 10^{-4}$, with an input measurement signal-to-noise ratio (SNR) of $20$ dB to ensure the robustness of the proposed models against noise in the testing data [34]. The proposed networks, DuRIN-$1$ and DuRIN-$2$, are trained on $1$-D data, and operate on $1$-D, $2$-D, and $3$-D data trace-by-trace. In the following sections, we provide experimental results on synthetic, simulated, and real data.

The Marmousi$2$ model [69] (width $\times$ depth: $17$ km $\times~{}3.5$ km) is an expansion of the original Marmousi model [74] (width $\times$ depth: $9.2$ km $\times~{}3$ km depth). The model provides a simulated $2$-D collection of seismic traces, and is widely used in reflection seismic processing for the calibration of techniques in structurally complex settings. The model has a $2$ ms sampling interval, with traces at an interval of $12.5$ m. We obtained the ground-truth reflectivity profiles $({\boldsymbol{x}})$ from the P-wave velocity and density models (1), and convolved them with a $30$ Hz Ricker wavelet to generate the input data $({\boldsymbol{y}})$, with measurement SNR of $20$ dB to test the robustness of the proposed networks to noisy conditions.

Here, we present results from a region of the Marmousi$2$ model that contains a gas-charged sand channel [69]. Table VII shows a comparison of the performance of the proposed DuRIN-$1$ and DuRIN-$2$ with that of the benchmark techniques for the portion containing the gas-charged sand channel. The Marmousi$2$ model [69] has a large number of low-amplitude reflections. Additionally, in Figure 5, we observe that BPI and SBL-EM introduce spurious reflections. We attribute the low PES scores observed for BPI and SBL-EM in Table VII to such spurious supports complementing the low-amplitude reflections in the Marmousi$2$ model. We re-evaluate the objective metrics after muting reflections with amplitudes $<~{}1\%$ of the absolute of the maximum amplitude, the results for which are reported in Table VIII and Figure 11. The re-evaluation did not affect the CC, RRE, and SRER scores adversely, but, as the PES is computed over significant reflections, it is higher (worse) for BPI and SBL-EM, and lower (better) for DuRIN-$1$ and DuRIN-$2$. The objective metrics reported in Table VIII show superior amplitude and support recovery accuracy of the proposed networks over the benchmark methods. Figure 11 shows that DuRIN-$1$ and DuRIN-$2$, along with SBL-EM, preserve the lateral continuity of interfaces, highlighted in the insets at the edge of the sand channel. The insets also show that BPI and FISTA introduce spurious false interfaces due to the interference of multiple events. DuRIN-$1$ and DuRIN-$2$ show limited recovery of the low-amplitude reflection coefficients right below the gas-charged sand channel (Figure 11 (f) and (g)), which could be investigated in future work.

To validate the proposed networks on real data, we use a $3$-D seismic volume from the Penobscot $3$D survey off the coast of Nova Scotia, Canada [70]. We present results from a smaller region of the $3$-D
volume, with $201$ Inlines ($1150$-$1350$) and $121$ Xlines ($1040$-$1160$), a region that includes the two wells of this survey (wells L-$30$ and B-$41$, [75]). Along the time axis, the region contains $800$ samples between the time interval of $0$ to $3196$ ms, with a $4$ ms sampling interval, and a $25$ Hz Ricker wavelet fits the data well [75].

Figure 12 shows the observed seismic data and recovered reflectivity profiles for Xline $1155$ of the survey, and overlaid in black, the seismic and reflectivity profiles, respectively, computed from the sonic logs of well L-$30$ [75]. The inverted reflectivity profiles for BPI and FISTA are smooth and lack detail. On the other hand, the predicted reflectivity profiles generated using SBL-EM, DuRIN-$1$ and DuRIN-$2$ provide more details for characterizing the subsurface. Additionally, DuRIN-$1$ and DuRIN-$2$ also resolve closely-spaced interfaces better and preserve the lateral continuity, evident from the interfaces around $1.1$ s on the time axis (Figure 12 (e) and (f)). These results demonstrate the capability of the proposed networks, namely, DuRIN-$1$ and DuRIN-$2$, on real data from the field. We note that the models are trained on synthesized seismic traces before being tested on the real data.