Seismic interpolation

Reconstruction of a 2D seismic gather from a random subset of its traces.

• True model: the fully sampled gather.

• Operator: pylops.Restriction along the spatial axis, keeping perc_subsampling of the traces (selected with a fixed seed).

• Data: the restricted gather, plus optional white Gaussian noise (noise_std, also drawn with a fixed seed).

Solvers compared:

Solver	Formulation
LSQR-SecondDerivative	Regularized least squares with a second derivative along the spatial axis, solved with LSQR
FISTA-FFT2D	L1 sparsity-promoting inversion in the f-k domain (2D FFT as preconditioner), solved with FISTA with a decaying sparsity weight
FISTA-Radon2D	L1 sparsity-promoting inversion in the linear Radon (tau-p) domain, solved with FISTA with a decaying sparsity weight
OMP-Radon2D	Greedy L0 inversion in the linear Radon domain, solved with Orthogonal Matching Pursuit
AIHT-Radon2D	L0-regularized inversion in the linear Radon domain, solved with Accelerated Iterative Hard Thresholding (PyProximal) with a decaying sparsity weight
POCS-L0-Radon2D	L0-regularized inversion in the linear Radon domain with exact data consistency, alternating hard thresholding with a decaying sparsity weight and projection onto the data-consistent affine set (PyProximal HQS)
OMP-FFT2D	Greedy L0 inversion in the f-k domain, solved with Orthogonal Matching Pursuit
AIHT-FFT2D	L0-regularized inversion in the f-k domain, solved with Accelerated Iterative Hard Thresholding (PyProximal) with a decaying sparsity weight
POCS-L0-FFT2D	L0-regularized inversion in the f-k domain with exact data consistency, alternating hard thresholding with a decaying sparsity weight and projection onto the data-consistent affine set (PyProximal HQS)

All Radon solvers use the same dictionary, pylops.signalprocessing.FourierRadon2D with a linear kernel. All f-k solvers use pylops.signalprocessing.FFT2D (adjoint) as the dictionary, with the restriction operator rebuilt in complex dtype and the real part taken after it.

Source: benchmarks/seismic_interpolation.

Results

Raw results: results/seismic_interpolation.csv.

Linear events

2D gather with three linear events, randomly subsampled in space.

The gather has 140 traces (2 m spacing) and 200 time samples (4 ms sampling). It contains three linear events with intercept times 0.2, 0.4 and 0.5 s, angles 0, 70 and -80 degrees (velocity 800 m/s) and amplitudes 1, -2 and 0.5, all with a 20 Hz Ricker wavelet.

Dataset parameters: noise_std=0, perc_subsampling=0.4, seed=42

Rank	Solver	MAE	PSNR [dB]	Time [s]	Figure
1	AIHT-FFT2D	0.002455	60.30	12.700	view
2	POCS-L0-FFT2D	0.002462	56.91	20.665	view
3	FISTA-FFT2D	0.002373	56.31	18.507	view
4	OMP-FFT2D	0.00454	54.61	23.115	view
5	AIHT-Radon2D	0.004103	52.86	43.547	view
6	FISTA-Radon2D	0.003084	52.32	33.120	view
7	POCS-L0-Radon2D	0.008483	50.12	161.054	view
8	OMP-Radon2D	0.005183	49.25	84.902	view
9	LSQR-SecondDerivative	0.007563	38.11	0.069	view