Theory and practice of interpolation in the pyramid domain

Putting everything together

We are now ready to present the missing data interpolation algorithm. One side effect of the fine sampling of the $u$ axis is that empty bins will appear in the pyramid domain. Missing data in ${\bf d}$ will add even more empty locations. We propose interpolating the missing data in ${\bf d}$ by filling the empty bins in ${\bf m}$. For this, we follow the approach of Claerbout and Fomel (2002). First, we want to honor the data where they are known by introducing the residual vector

$${\bf r_d} = {\bf W_d(Lm - d)},$$ (13)

where ${\bf W_d}$ is a masking operator equal to unity where data are known, and zero where data are missing. Solving for ${\bf m}$ minimizing the amplitude of ${\bf r_d}$ only will not fill the empty bins. We need to add a regularization term that will enforce a certain multivariate spectrum to the vector ${\bf m}$:

$${\bf r_m} = {\bf Am},$$ (14)

where ${\bf A}$ is a pef in the model space. Assuming that the pef ${\bf A}$ is known, we can fill the empty locations in ${\bf m}$ by minimizing $f({\bf m})=\Vert{\bf r_d}\Vert _2+\epsilon \Vert{\bf r_m}\Vert _2$, where $\epsilon$ is a balancing operator between data fitting and model space regularization.

There are two issues with this approach. The first issue is that the convergence towards a solution will be slow. To accelerate this process we introduce a new variable ${\bf q}={\bf Am}$ and rewrite equations (13) and (14) as follows:

$$\begin{array}{lcl} {\bf r_d} & = & {\bf W_d(LA^{-1}q - d)} \\ {\bf r_q} & = & {\bf q} \end{array}$$ (15)

We then minimize $f({\bf q})=\Vert{\bf r_d}\Vert _2+\epsilon \Vert{\bf r_q}\Vert _2$ and compute ${\bf\tilde m}={\bf A^{-1} \tilde q}$, where ${\bf\tilde q}$ minimizes $f({\bf q})$. The term ${\bf A^{-1}}$ is computed by applying a polynomial division to ${\bf q}$ which yields fast filling of the empty bins. Because ${\bf A}$ is a miminum-phase filter, the polynomial division is stable and will not cause the solution to blow-up. This preconditioning of the problem has been used in numerous geophysical problems [Herrmann et al. (2009); Fomel and Guitton (2006); Clapp et al. (2004)]. In practice, we set $\epsilon=0$ and minimize ${\bf r_d}$ only [Trad et al. (2003); Guitton and Claerbout (2004)].

The second issue is that both ${\bf q}$ and ${\bf A}$ are unknown in equation (15). To circumvent this problem we bootstrap the pef estimation by first assuming that ${\bf A}$ is a 1-D gradient. To make sure that we can apply the polynomial division, we set the second coefficient of the gradient to $-0.96$ instead of $-1$. We then minimize $f({\bf q})$ with this first pef, find a new ${\bf\tilde m}$ and estimate a better pef ${\bf A}$ from it. Having a better pef, we can minimize $f({\bf q})$ again:

		 ${\bf A} \leftarrow (1 -0.96)^{T}$ 


iterate {                   


minimize $f({\bf q}$)       

${\bf\tilde m} \leftarrow {\bf A^{-1} \tilde q}$ 


estimate ${\bf A}$ from ${\bf\tilde m}$

}

We are essentially solving the non-linear problem in a step-wise fashion by keeping ${\bf A}$ constant within each non-linear loop. In practice, we notice that only 4 to 5 non-linear iterations are necessary to converge towards a pef that yields accurate interpolation of the missing data.

In the next section, we apply this algorithm to synthetic and field data examples in 2-D. We show that aliased and irregularly-sampled data can be interpolated.

Theory and practice of interpolation in the pyramid domain