Seismic data interpolation using streaming prediction filter in the frequency domain

$f$ -$x$ streaming prediction filter

In the $f$ -$x$ domain, the PF predicts the data along the spatial direction, and the relationship between data and prediction filter can be summarized as

$$\begin{aligned}g_{m,n} & \approx \sum_{p} g_{m,n-p} f_{m,n,p} \... ...rox \mathbf{G}_{m,n}^{\mathsf{T}} \mathbf{F}_{m,n}, \end{aligned}$$

where $m$ and $n$ are the indices of the seismic data sample $g_{m,n}$ along the frequency $f$ axis and space $x$ axis, respectively. Vector $\mathbf{F}_{m,n}$ is the group of filter coefficients $f_{m,n,p}$ in the adaptive prediction filter, and each group $\mathbf{F}_{m,n}$ corresponds to a data sample $g_{m,n}$ . Vector $\mathbf{G}_{m,n}$ denotes several data points $g_{m,n-p}$ with spatial shift $p$ near $g_{m,n}$ . Spatial shift $p$ is related to the filter size (the number of filter coefficients in $\mathbf{F}_{m,n}$ ), and the filter size should theoretically be larger than or equal to the number of seismic events contained in the local space window. When the causal filter structure is considered, as shown in Fig. 1a, the spatial shift is chosen as $p \in [1, 3]$ , vector $\mathbf{G}_{m,n}$ can be expressed as $\mathbf{G}_{m,n} = \{ g_{m,n-1}, g_{m,n-2}, g_{m,n-3} \}$ , and $\mathbf{F}_{m,n} = \{ f_{m,n,1}, f_{m,n,2}, f_{m,n,3} \}$ . For the non-causal filter structure (Fig. 1b), the spatial shift is $p \in [-3,-1] \cup [1, 3]$ , and vector $\mathbf{G}_{m,n}$ and $\mathbf{F}_{m,n}$ can be represented as $\mathbf{G}_{m,n} = \{ g_{m,n+3}, g_{m,n+2}, g_{m,n+1}, g_{m,n-1}, g_{m,n-2}, g_{m,n-3} \}$ , $\mathbf{F}_{m,n} = \{ f_{m,n,-3}, f_{m,n,-2}, f_{m,n,-1}, f_{m,n,1},\ f_{m,n,2}, f_{m,n,3} \}$ .

With Eq. (3), we established the following minimization problem, like Eq. (1), to calculate filter $\mathbf{F}_{m,n}$ :

$$\min_{\mathbf{F}_{m,n}} \Vert g_{m,n} - \mathbf{G}_{m,n}^{\mathsf{T}} \mathbf{F}_{m,n} \Vert _{2}^{2}.$$ (4)

According to the above explanation of vector $\mathbf{F}_{m,n}$ , there are several unknown filter coefficients, yet we established only one equation. Eq. (4) describes an ill-posed problem, which requires constraints to obtain a stable solution. The framework of the streaming computation (Fomel and Claerbout, 2016; Sacchi and Naghizadeh, 2009; Liu and Li, 2018) establishes the constraint relationship by using local smoothness. We extended this method into the frequency domain, including the $f$ -$x$ domain and the $f$ -$x$ -$y$ domain (the next section) to stabilize the filter coefficient solution. Additionally, we discussed in detail the effects of filter structure and processing path, and provided the corresponding interpolation algorithms. Here, multiple constraints based on local smoothness are used to constrain the solution of Eq. (4):

$$\begin{aligned}\min_{\mathbf{F}_{m,n}} \Vert g_{m,n} & - \mathb... ...athbf{F}_{m,n} - \mathbf{F}_{m,n-1} \Vert _{2}^{2}, \end{aligned}$$

where $\lambda_{f}$ and $\lambda_{x}$ are the weights of the regularization terms along the frequency $f$ and space $x$ axes. $\lambda_{f}^{2} \Vert \mathbf{F}_{m,n} - \mathbf{F}_{m-1,n} \Vert _{2}^{2}$ shows the local smoothness along the frequency direction ( $\lambda_{f}\mathbf{F}_{m,n} \approx \lambda_{f}\mathbf{F}_{m-1,n}$ ). Likewise, $\lambda_{x}^{2} \Vert \mathbf{F}_{m,n} - \mathbf{F}_{m,n-1} \Vert _{2}^{2}$ controls the local smoothness along the space direction as $\lambda_{x}\mathbf{F}_{m,n} \approx \lambda_{x}\mathbf{F}_{m,n-1}$ . The block matrix Eq. (6) has the same solution as Eq. (5), and it demonstrates the effect of the local smoothness constraints; when $\mathbf{F}_{m-1,n}$ and $\mathbf{F}_{m,n-1}$ are considered known, the local smoothness conditions ( $\lambda_{f}\mathbf{F}_{m,n} \approx \lambda_{f}\mathbf{F}_{m-1,n}$ and $\lambda_{x}\mathbf{F}_{m,n} \approx \lambda_{x}\mathbf{F}_{m,n-1}$ ), as newly added equations, can be used to stabilize the solution of $\mathbf{F}_{m,n}$ :

$$\begin{bmatrix}\mathbf{G}_{m,n}^{\mathsf{T}} \ \lambda_{f} \math... ...\lambda_{f} \mathbf{F}_{m-1,n} \ \lambda_{x} \mathbf{F}_{m,n-1} \end{bmatrix}.$$ (6)

One can obtain the following least-squares solution of Eq. (5) and (6):

$$\begin{aligned}\mathbf{F}_{m,n}= & {( \lambda_{f}^{2} \mathbf{I... ...f{F}_{m-1,n} + \lambda_{x}^{2}\mathbf{F}_{m,n-1} ), \end{aligned}$$

where $\{\bullet\}^{*}$ denotes the conjugate operator. Let

$$\begin{cases}\lambda^{2} = \lambda_{f}^{2} + \lambda_{x}^{2} ... ...hbf{F}_{m-1,n} + \lambda_{x}^{2}\mathbf{F}_{m,n-1} \end{cases},$$ (8)

we get a simplified equation:

$$\mathbf{F}_{m,n} = {( \lambda^{2} \mathbf{I} + \mathbf{G}_{m,n}^{... ...}^{-1} ( g_{m,n} \mathbf{G}_{m,n}^{*} + \lambda^{2} \mathbf{\tilde{F}}_{m,n} ).$$ (9)

Meanwhile, $(\lambda^{2}\mathbf{I} + \mathbf{G}_{m,n}^{*}\mathbf{G}_{m,n}^{\mathsf{T}})$ has the analytical inversion $\frac{1}{\lambda^{2}}(\mathbf{I} - \frac{\mathbf{G}_{m,n}^{*}\mathbf{G}_{m,n}^... ...hsf{T}}} {\lambda^{2} + \mathbf{G}_{m,n}^{\mathsf{T}}\mathbf{G}_{m,n}^{*} } )$ when one extends the Sherman-Morrison formula (Sherman and Morrison, 1950; Hager, 1989; Bartlett, 1951) to the complex space. We can obtain the analytical solution of the $f$ -$x$ SPF:

$$\begin{aligned}\mathbf{F}_{m,n} & = \frac{1}{\lambda^{2}}(\math... ...hsf{T}}\mathbf{G}_{m,n}^{*} } \mathbf{G}_{m,n}^{*}. \end{aligned}$$

Eq. (10) shows a recursive relationship from the previous filters ( $\mathbf{F}_{m-1,n}$ and $\mathbf{F}_{m,n-1}$ ) to current filter $\mathbf{F}_{m,n}$ . For data interpolation, Eq. (2) becomes a well-posed problem, and the unknown data sample can be calculated by

$$\hat{g}_{m,n} = \mathbf{G}_{m,n}^{\mathsf{T}} \mathbf{F}_{m,n}.$$ (11)

In 2D data interpolation, we proposed the following Algorithm 1 to reconstruct the missing seismic data by using the 2D $f$ -$x$ SPF. To start, $\mathbf{F}_{m-1,n}$ and $\mathbf{F}_{m,n-1}$ are initialized to $\mathbf{0}$ . When the processing path in Algorithm 1 is followed, both $\mathbf{F}_{m-1,n}$ and $\mathbf{F}_{m,n-1}$ are known. We designed a space-causal filter form (Fig. 1a) for the $f$ -$x$ SPF. Fig. 1b suggests that the space-noncausal form may involve the interference of the unknown data samples, but the causal one can avoid this problem.

causal2d,noncausal2d
Figure 1. Space-causal filter form (a) and space-noncausal filter form (b) in the $f$ -$x$ domain. Solid circle denotes known sample, and hollow circle denotes unknown sample.

Seismic data interpolation using streaming prediction filter in the frequency domain