Adaptive subtraction

SRME produces a multiple prediction that is close but not exactly right: the amplitude may be 5 to 30 % off, the timing may be off by a few milliseconds, and the wavelet shape is not quite the true multiple wavelet. Direct subtraction would leave a ringing remnant at every multiple location, which is why SRME always ends with an adaptive subtraction: fit a short filter to absorb the mismatch before subtracting.

1. The problem direct subtraction leaves

Observed: $d = p + m$ (primary plus true multiple). Predicted: $\\hat m \\approx m$, but wrong in amplitude, timing and shape. Direct residual:

The $m - \\hat m$ term is not small: with a timing error of a few milliseconds it can carry more energy than the multiple itself, and on a stack it appears as a ringing copy of every multiple. In the figure, change how the prediction is wrong, then design the filter that corrects it and watch what each subtraction leaves:

The prediction opens at 0.79 of the true amplitude (0.55 against 0.70) and 15 ms late, about half a lobe of the 28 Hz wavelet. Subtracted as it stands it leaves 4.1 dB more energy than the multiple had, which is what the closed form $1 + a^2 - 2a\\rho(\\tau)$ predicts, with $\\rho(\\tau)$ the normalised correlation of the wavelet with itself at the time error. A 21-point two-sided matching filter takes the multiple down by 17.4 dB, and its peak, refined between the coefficients at −16 and −12 ms, sits at −15.1 ms, the advance the time error needs. The ceiling is the 3 % noise in the design window: the filter fits some of it too, so even an exact prediction, once matched, stops near 17 dB, and 11 points (19.6 dB) do better than 61 (15.4 dB).

2. The Wiener filter formulation

Find a short two-sided filter $f$ of length $L$ that minimizes $E = \\sum\_{t \\in W} \\left(d\_t - (f \* \\hat m)\_t\\right)^2$ over a design window $W$. Setting $\\partial E/\\partial f = 0$ yields

where $R\_{\\hat m\\hat m}$ is the Toeplitz autocorrelation of the predicted multiple in the design window, and $R\_{d\\hat m}$ is the cross-correlation of the observed data with the prediction, at lags on both sides of zero. Solve it by Levinson recursion, as the figure does: it is the algebra of spiking decon (Section 2.6) with a different goal. The two-sided lags matter. A causal filter can only delay, so against the late prediction in the figure it manages only 4.6 dB.

Apply $f$ outside the design window too: compute $f \* \\hat m$ over the whole trace and subtract it from $d$. This works where the mismatch is stationary. Where it is not, use time-varying windows (see Advanced variants below).

3. Tuning the filter

• Length. Short (5 to 10 samples) absorbs amplitude errors and small time shifts. Long (30 to 60 samples) captures richer wavelet-shape mismatches. Too long, and it starts fitting noise and primary energy too: over-subtraction. With a primary under the multiple in the figure, 11 points take 42 % of its amplitude and 61 points take 87 %, yet even 3 points take 19 % while removing under 2 dB of the multiple: shorter filters take less of the primary and less of the multiple, and where a primary overlaps the multiple no length is safe.
• Design window. It must be centred on a region where the multiple dominates and the primary is weak. If primaries overlap the multiple inside the window, the filter will fit part of them and subtract it, a long filter more than a short one. A window that holds no predicted multiple gives the filter nothing to fit at all.
• Window width. Narrow catches local phase variation; wide averages. 100-300 ms is typical.
• White-noise stabilization. Same trick as decon: add $\\varepsilon\\, r\_0$ (with $\\varepsilon$ about 0.1 to 1 %) to the diagonal of the autocorrelation so the solve stays stable.

4. The production recipe

1. Predict the multiple (via SRME, predictive decon, or Radon).
2. Define design windows where the multiple dominates (picked manually or by time gates below strong reflectors).
3. Solve for $f$ per window, per trace (or per gather, with a single filter shared).
4. Convolve $f \* \\hat m$ and subtract it from the data.
5. QC: check that the primary is preserved and the multiple is attenuated. Repeat with different windows if over-subtraction occurred.

5. Advanced variants

• Time-varying adaptive. A different filter for every few hundred milliseconds of the trace handles a non-stationary wavelet.
• Pattern-based (PEF) adaptive. Instead of a least-squares (Wiener) matched filter, use prediction-error filters to capture the local pattern of the multiples and subtract them adaptively. Better where multiples and primaries overlap in time.
• Transform-domain adaptive. Separate primaries and multiples in a sparsifying domain (curvelet or wavelet), where they carry distinct signatures. Robust to dispersion and dispersive phase errors in the prediction.
• L1-regularized adaptive. Sparsity-promoting filter design; less prone to over-subtraction in areas where the prediction is weak.
• Multi-trace adaptive. Fit a single filter to all traces in a CMP gather; enforces spatial smoothness of the correction.

6. The one failure mode to watch for

Over-subtraction. If a primary overlaps the predicted multiple inside the design window, the filter can reshape the prediction to fit part of it, and subtracting that removes part of the primary, leaving it dimmed or holed. A primary alone in a window with no predicted multiple is safe: there is nothing for the filter to fit. Guard against the overlap with short filters, windows placed where the multiple dominates, and a QC of primary amplitudes before and after (a difference plot).

Where this goes next

Section 4.5 closes Part 4 with the hardest multiple class: inter-bed multiples that never reach the surface. SRME cannot predict them because they never reach the surface, and they often have too little moveout for Radon. They are predicted either from a model of the overburden or with data-driven methods such as the inverse scattering series, and then removed with adaptive subtraction again.

References

• Verschuur, D. J., Berkhout, A. J., Wapenaar, C. P. A. (1992). Adaptive surface-related multiple elimination. Geophysics, 57, 1166.
• Berkhout, A. J., Verschuur, D. J. (1997). Estimation of multiple scattering by iterative inversion. Geophysics, 62, 1586.
• Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
• Claerbout, J. F. (1976). Fundamentals of Geophysical Data Processing. McGraw-Hill.