Modeling data error during deconvolution

INTRODUCING NOISE AS ITSELF A MODEL

The idea of this paper is that we should not try push all our data into the convolutional model. We should explicitly solve for an unknown part of the data that poorly fits this model. I call this part noise and define it negatively $-N$ (so the minus sign is missing from all the analysis and code).

The decon filter $C=e^U$ , parameterized by $U$ , we take as noncausal. The constraint is no longer a spike at zero lag, but a filter whose log spectrum vanishes at zero lag, $0 = u_0 = \int \ln C(\omega) \ d\omega$ , so we are now constraining the mean of the log spectrum. This is a fundamental change which we confess to being somewhat mysterious.

The single regression for $U$ including noise $N$ now becomes two.

$$\approx_h (D+N)e^U = (D+N)C$$ 0 (1)

$$\approx_2 N$$ 0 (2)

The notation $\approx_h$ means the data fitting is done under a hyperbolic penalty function. The regularization need not be $\ell_2$ . To save clutter I leave it as $\ell_2$ until the last step when I remind how it can easily be made hyperbolic.

Under the constraint of a causal filter with $c_0=1$ , traditional auto regression for $c_t= {\rm FT}^{-1}C$ with its regularization looks more like

$$\approx N\ =\ Dc$$ 0 (3)

$$\approx c$$ 0 (4)

Comparing equations 1-2 with 3-4 you see we are not simply rehashing traditional methodology but seem to be off in a wholly new direction! We are here because $C=e^U$ solved our non-minimum phase problem, and seeing sea swell in our estimated shot wavelets told us we need to replace $D$ by $D+N$ .

Antoine noticed the quasi-Newton method of data fitting requires gradients but not knowledge of how to update residuals $\Delta \bold r$ so the only thing we really need to think about is getting the gradient. The gradient wrt $U$ is the same as before (Claerbout et al. (2011)) except that $D+N$ replaces $D$ . The gradient wrt $N$ is the new element here.

Let $d$ , $n$ , and $c$ be time functions (data, noise, and filter). Let $r=(d-n)\ast c$ be the residual. Let $h_t=h(r_t) = {\rm hyperbolic\ stretch\ of\ } r$ . Expressing our two regressions in the time domain we minimize

$$\min_n \quad \sum_t n^2/2 \ +\ h((d+n)\ast c)$$ (5)

A scaling factor is required between the terms. We expect to learn it by experimentation.

Now we go after the gradient, the derivative of the penalty function wrt each component of noise $n_s$ . Let the derivative of the penalty function $h(r_t)$ wrt its argument $r_t$ be called the softclip and be denoted $h'_t = h'(r_t)$ . Let $H'$ denote the FT of $h'$ . Let $c'(t)$ be the time reverse of $c(t)$ while in Fourier space $C'$ is the conjugate of $C$ .

$$\Delta n_s = n_s \ +\ \ \frac{\partial}{\partial n_s} \ \sum_t \ h(r_t)$$ (6)

$$= n_s \ +\ \ \sum_t h'(r_t) \ \frac{\partial}{\partial n_s} \ (d+n) \ast c$$ (7)

$$= n_s \ +\ \ \sum_t h'_t \ \frac{\partial}{\partial n_s} \ \sum_\tau n_{t-\tau} c_\tau$$ (8)

$$= n_s \ +\ \ \sum_t h'_t \ \frac{\partial}{\partial n_s} \ \sum_j n_{j} c_{t-j}$$ (9)

$$= n_s \ +\ \ \sum_t h'_t \ c_{t-s}$$ (10)

$$= n_s \ +\ \ \sum_t h'_t \ c'_{s-t}$$ (11)

$$\Delta N = N \ +\ C' H'$$ (12)

For simplicity I set out with a quadratic penalty function on the noise, but it is easy to make it hyperbolic. Simply use softclip on $n$ . Change $\Delta n = n+\cdots$ to $\Delta n = h'(n)+\cdots$ .

Now having the gradient we should be ready to code.

Modeling data error during deconvolution