Data-driven internal multiple elimination by the Marchenko method (2D)

Theory

Multiple Modeling Marchenko 2D removes internal (interbed) multiples from the shot records of a 2D line with the Marchenko multiple elimination (MME) method. The multiples are predicted from the recorded reflection response itself: no velocity model, no picked horizons and no adaptive subtraction are needed, and all orders of internal multiples from all generators are handled in one process.

For every selected shot the module solves the projected Marchenko equations as a Neumann series - an iterative sum of repeated correlations and convolutions with the reflection response R - for each time level inside a truncation window. The result is the shot record with the internal multiples removed (primaries). The multiples are the difference between the input record and the primaries.

Two schemes are available:

• MME - internal multiple elimination; the primaries keep their original amplitudes.
• T-MME - additionally compensates the primaries for transmission losses in the overburden.

The module works in the frequency domain on a regular grid along the line. The source and receiver positions are projected onto a best-fit straight line and placed on a common grid with half the station spacing. Both half-station staggered acquisition (sources and receivers on alternate grid slots) and co-located acquisition (sources at receiver positions) are supported; on staggered geometry the data are resampled between the two station families with a 16-tap interpolator.

This is a batch module: it reads R from disk and writes the results to new data files. The whole reflection response is held in memory during the run, so a real line can need tens of gigabytes of memory and several hours of computing time. Use Analyze geometry and a test run on a few shots before starting a full run.

Input data

Reflection response R

Prestack shot records of the 2D line (data handle and trace headers). R is the only seismic input: the shots to demultiple are selected from it, and the output is written in its amplitude domain. Source and receiver coordinates must be filled in the trace headers; the SOURCE_SP header identifies the shots.

R must be prepared before it is used. It is neither a stack nor a raw shot gather:

• the direct wave must be removed;
• free-surface (surface-related) multiples must be eliminated, e.g. with SRME;
• the source wavelet must be deconvolved (designature);
• the amplitudes must be calibrated so that R represents dimensionless reflection coefficients of order 1.

If raw or incorrectly scaled data are used, the method does not necessarily fail: it can converge to a plausible-looking but wrong result, or diverge (see Information).

Source wavelet

Wavelet gather; the first trace is used, and its centre sample is taken as time zero. The wavelet is convolved into the initial estimate only, to give the output the bandwidth of the data. It does not deconvolve R.

Output primaries file

Data file (.gsd) for the demultipled shot records (primaries).

Output multiples file

Data file (.gsd) for the predicted internal multiples (input record minus primaries). At least one of the two output files must be set. Both outputs contain the selected shots of R with their original trace headers and sample interval.

Parameters: Shot selection

First shot (SOURCE_SP)

SOURCE_SP value of the first shot to process. Default: 0.

Last shot (SOURCE_SP)

SOURCE_SP value of the last shot to process. Default: 0. When both First shot and Last shot are 0, all shots of R are processed. A small range is recommended for test runs.

Shot decimation

Process every N-th selected shot (minimum 1, default 1 = every shot). N counts shots, not SOURCE_SP units, so gaps in the SOURCE_SP numbering do not cause uneven sampling. Selected shots whose source position has no matching source in R are skipped and reported in the log.

Parameters: Computation

Min frequency

Lowest frequency used in the computation, in Hz. Default: 0 Hz.

Max frequency

Highest frequency used, in Hz. Default: 70 Hz. Must be greater than Min frequency. Set the band to the useful bandwidth of the data; a narrower band also reduces memory use and run time.

Grid step (0 = auto)

Spacing of the common grid along the line, in metres. Default: 0, meaning automatic (half of the station spacing). Enter a value only if Analyze geometry shows a wrong step, slot collisions or a BROKEN acquisition.

Lateral taper

Options: Off, Shot, R, Both. Default: Off. In the current version this parameter has no effect.

Taper stations

Default: 0. In the current version this parameter has no effect.

Parameters: Marchenko iterations

Number of iterations

Number of Neumann-series terms computed for each time level (minimum 1, default 10). More iterations predict higher-order multiples more completely but take longer. At least 2 are needed for the convergence check described under Information.

Start time

First time level computed, in seconds. Default: 0.16 s. Only times between Start time and End time are processed: outside this range the primaries output contains no energy, and the multiples output equals the input.

End time (0 = trace end)

Last time level computed, in seconds. Default: 0 = end of the trace. Must be greater than Start time when set.

Transmission compensation (T-MME)

Off (default): MME scheme. On: T-MME scheme, which also compensates the primaries for transmission losses. The selected scheme is written to the log.

Parameters: Mute window

Shift (wavelet length)

Length of the wavelet, in samples (default 20). It sets the truncation window of the Marchenko equations: the window starts Shift samples after zero time, and its end follows the current time level (Shift samples before it for MME, Shift samples after it for T-MME). The window is empty for time levels below 2 × Shift samples. Set it to the length of the wavelet in R.

Smooth

Length of the cosine taper on the edges of the truncation window, in samples. Default: 10.

Parameters: Reflection response correction

These corrections are applied to every trace of R as it is read.

Scale

Scale factor of the reflection response. Default: 2.0. The factor 2 is required by the projected Marchenko equations - do not change it. If the amplitudes of R are wrong, apply a constant gain to the input data instead.

t^n gain exponent

Multiplies R by tⁿ (t = time in seconds) to compensate geometric spreading if this has not been done in preprocessing. Default: 0 (off).

Q

Quality factor for constant-Q absorption compensation, which multiplies R by exp(π · f0 · t / Q). Default: 0 (off). Q and Q reference frequency must both be set, or both left at 0. A small Q with a high reference frequency over a long record can overflow; this is detected and reported.

Q reference frequency

Reference frequency f0 of the Q compensation, in Hz. Default: 0 (off).

Pad

Number of zero samples added to the traces before the Fourier transform, to avoid wrap-around. Default: 1024.

Information

Geometry summary

One-line summary of the geometry used by the run: grid step, number of grid slots, number of source and receiver stations, slot collisions and the acquisition type. It is updated by Analyze geometry and at the start of every run. The full analysis is written to the log.

The log also records the full parameter set at the start of every run, and, after the first shot, the ratio of the predicted result to the input and a convergence probe: the peak operator norm over the computed time levels, which must be below 1.

The run is stopped with an error message in these cases:

• Stations do not sit consistently on alternating grid slots (BROKEN acquisition) - the grid step is wrong for this line. Run Analyze geometry and check the reported step and residual.
• Every selected shot has fewer than 16 live receiver stations - the output would be interpolated from zeros. Select shots with a longer aperture. If only some shots are short (usually at the ends of the line), a warning is written and their multiple model should not be trusted.
• The initial guess is all zero - the first shot has no live traces in R.
• The multiple prediction is identically zero - sources and receivers are not related on the common grid; check Analyze geometry.
• The Neumann series diverged - usually the amplitudes of R are too large. The message gives the measured operator norm and, when it is above 1, a suggested constant gain for the input data. Apply that gain to R and run again; do not change Scale.

Actions

Analyze geometry

Fits the line and builds the common grid without running the solver; no output files are needed. The log lists the number of traces and stations, the cross-line deviation of the stations from the fitted line (rms and maximum), the grid step (and the automatic value), the number of grid slots, the largest station-to-slot residual, source and receiver slot collisions, the slot parity of both station families and the acquisition type: half-station staggered, co-located or BROKEN. A BROKEN geometry will be refused by the run. Always run this action first, because a full run takes hours.

Examples

Recommended workflow:

• Prepare R: remove the direct wave, eliminate surface-related multiples, deconvolve the source wavelet and calibrate the amplitudes.
• Connect R and the source wavelet, and run Analyze geometry. Check that the acquisition is staggered or co-located and that there are no slot collisions.
• Make a test run on a few shots (First shot / Last shot, or Shot decimation). Check the convergence probe in the log and compare the primaries and multiples outputs with the input.
• Run the whole line and QC the results.

Related articles

• IMP (Interbed multiple prediction)
• Interbed multiple modeling
• SRME 2D/3D

References

Thorbecke, J., Zhang, L., Wapenaar, K., and Slob, E., 2021, Implementation of the Marchenko multiple elimination algorithm: Geophysics, 86(2).

* * *   If you have any questions, please send an e-mail to: support@geomage.com  * * *