Spiking deconvolution (Wiener filtering)

The convolutional model of Section 0.3 says $s = w \\ast r + n$: the trace is the wavelet convolved with the reflectivity, plus noise. If we knew the wavelet, we could invert it to recover $r$. The catch: we do not. Spiking deconvolution estimates an operator that approximately undoes the wavelet, using only the trace itself.

1. What “spiking” means

Design a filter f\[n\] such that

that is, the operator convolved with the wavelet produces a spike. Apply $f$ to the trace: $f \\ast s = f \\ast w \\ast r + f \\ast n \\approx \\delta \\ast r = r$, as long as the filtered noise $f \\ast n$ stays small. You have the reflectivity.

The problem: $w$ is unknown. The trick: under the white-reflectivity assumption (reflection coefficients uncorrelated and zero-mean), the autocorrelation of the trace is proportional to the autocorrelation of the wavelet, in expectation and over a long enough window; white noise adds only to its zero lag. We compute that autocorrelation from the trace, then solve for $f$ with the Wiener-Hopf equations.

2. The Wiener-Hopf equation

The operator $f$ of $N$ samples (an operator length $L = N \\Delta t$) that minimizes $\\|f \\ast w - \\delta\\|\_2^2$ solves

where $R\_{ww}$ is the $N \\times N$ symmetric Toeplitz autocorrelation matrix and $e\_1 = (1, 0, \\dots, 0)^T$. (Strictly the right side is $w\_0 e\_1$; the unknown $w\_0$ only rescales $f$.) Under white reflectivity, replace $R\_{ww}$ with $R\_{ss}$, the matrix of the trace’s autocorrelation $\\phi\_k$. This is a Toeplitz solve: Levinson’s recursion does it in $O(N^2)$ operations instead of $O(N^3)$.

3. White-noise stabilization

The autocorrelation matrix is often nearly singular: the wavelet has a notch, or the signal is band-limited, and the inverse then asks for enormous gain where the spectrum is nearly empty. The standard fix is the white-noise percentage, or prewhitening: add a small fraction $\\varepsilon$ of the zero-lag autocorrelation to itself:

Adding $\\varepsilon \\phi\_0$ to the zero lag is the same as adding white noise of that power to the trace, so it puts a floor under the spectrum the operator inverts: the operator flattens the spectrum down to that floor and no further. Typical values are 0.1 % to 5 %. A higher $\\varepsilon$ gives a gentler operator, less compression and less noise; a lower one gives a sharper spike and lifts more noise. Noise already in the data does the same job: at a signal-to-noise ratio of 20 dB it adds about 1 % to the zero lag on its own, and below that level lowering $\\varepsilon$ changes the operator little: at 20 dB the boost of the operator stops at +28 dB, where on almost clean data it climbs to +34 dB.

In Figure 2.6 you design the operator yourself. Lower $\\varepsilon$ and see how far the wavelet collapses and what that costs in noise, shorten the operator, then break the two assumptions the method rests on: a minimum-phase wavelet and white geology.

At the starting $\\varepsilon$ of 10 % the operator widens the wavelet’s band from 36 Hz to 76 Hz. Lower it to 1 % and the band reaches 100 Hz: in (a) the ringing minimum-phase wavelet collapses to a short pulse, in (c) the output spectrum stays within about 6 dB of its peak from 5 Hz to 125 Hz, and in (d) the output lines up with the reflectors, its correlation with the true reflectivity rising from 0.35 for the recorded trace to 0.81. Spiking and whitening are one statement. The cost is noise: the operator lifts the weak high frequencies whether they carry signal or noise, and the output S/N falls from 30 dB to 16 dB. A longer operator and a smaller $\\varepsilon$ give sharper spikes and a higher noise floor at the band edges; a shorter operator and a larger $\\varepsilon$ are stable but leave the wavelet only partly compressed.

The figure also shows the method’s two failures. Switch the wavelet to zero phase: the trace looks different, but its autocorrelation is the same, so the operator is the same and the output autocorrelation is just as white. Yet $f \\ast w$ is no longer one-sided: 48 % of its energy arrives before time zero and only 23 % at it, against 66 % for the minimum-phase wavelet, and the output matches the true reflectivity at only 0.47, worse than the 0.70 of the recorded trace, because the operator assumed a minimum-phase wavelet. Switch the geology to thin beds and the operator reads their 28 ms top-and-base correlation as part of the wavelet: $f \\ast w$ in (a) carries an echo of half its peak at the bed period, which lands on each bed’s base with the opposite sign and cancels part of it. The pulse itself stays sharp (its first 16 ms have a band of 95 Hz, against 99 Hz over white geology), but the echo combs the spectrum, so the bandwidth readout falls from 100 Hz to 39 Hz. Shorten the operator to 20 ms, less than the bed period, and the echo all but goes.

4. The minimum-phase assumption

Wiener spiking decon produces a minimum-phase operator: all zeros of the inverse filter lie on the minimum-phase side of the unit circle (Section 0.6). If the true wavelet is also minimum phase, the operator works well. If the wavelet is zero phase (a Ricker, or vibroseis data after correlation) or mixed phase (most real data), spiking decon still compresses the wavelet but distorts its phase: events come out asymmetric, shift by a few milliseconds and can reverse polarity. The autocorrelation carries no phase, so nothing in the trace or in the output’s autocorrelation warns you.

Three tools come up in practice, and only one of them fixes the phase:

• Pre-whitening (what $\\varepsilon$ already does) stabilizes the inversion and caps how hard the operator boosts weak frequencies; it does not remove the phase error of a non-minimum-phase wavelet.
• Zero-phase post-processing applies a phase correction to restore zero-phase.
• Surface-consistent deconvolution (Section 2.8) decomposes the operator into source, receiver, offset and CMP terms, which averages out noise and gives consistent operators across the survey; it still assumes minimum phase.

5. Practical parameter picking

• Operator length. Long enough to capture the wavelet’s autocorrelation out to where it has essentially died (often 100 to 200 ms), and no longer than the design window can support: a very long operator starts fitting the geology and the noise in the autocorrelation estimate. Too short under-compresses.
• White noise $\\varepsilon$. Start at 1 %. Increase it if the output is noisy or rings, decrease it if the wavelet is not compressed enough, and remember that below the noise level of the data it does little.
• Design window. Estimate $R\_{ss}$ only over a time range where primaries dominate (no ground roll, no strong multiples). A design gate of about 2 to 3 s that covers the primaries of interest, below the first breaks and ground roll, is typical.
• Statics and decon. A pure time shift leaves a trace’s autocorrelation unchanged, so statics do not alter a trace-by-trace spiking operator. Field statics are often applied first so that the design gate sits in the same place on every trace; residual statics come after decon.

6. When spiking fails

• Non-white reflectivity. Cyclic reservoirs (thin-bed stratigraphy, rhythmic turbidites) violate the whiteness assumption. $R\_{ss} \\neq R\_{ww}$, so the operator is biased: it treats part of the geology as wavelet and removes it. Use a design gate that averages over a long window of geologically diverse sediments, or keep the operator well short of the cycle.
• Non-minimum-phase wavelet. Spiking compresses but distorts phase. Use zero-phase post-processing or surface-consistent decon (Section 2.8).
• Coloured noise inside the design window. Noise autocorrelation leaks into $R\_{ss}$ and the operator tries to whiten the noise as well. Mask noisy regions out of the design gate.
• Strong multiples inside the design window. Multiples have their own predictable autocorrelation, which spiking will partially invert. Use predictive decon (Section 2.7) explicitly.

Where this goes next

Section 2.7 uses the same Wiener framework for a different goal: predictive deconvolution. Instead of spiking the wavelet, we predict and subtract the part of the trace that repeats at a fixed lag, the multiple train from a water bottom or a strong reflector.

References

• Robinson, E. A. (1957). Predictive decomposition of seismic traces. Geophysics, 22, 767.
• Treitel, S., Robinson, E. A. (1966). The design of high-resolution digital filters. IEEE Trans. Geosci. Electron., 4, 25.
• Robinson, E. A., Treitel, S. (2008). Digital Imaging and Deconvolution. SEG.
• Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
• Claerbout, J. F. (1976). Fundamentals of Geophysical Data Processing. McGraw-Hill.