Encoded FWI & computational strategies

Production FWI is computationally extreme: a full-physics 3D acoustic FWI at 10 Hz on a modest survey (10 km × 10 km × 8 km, $\\Delta x = 25\\ \\text{m}$) takes thousands of shots, each requiring a forward plus an adjoint simulation per iteration, for tens of outer iterations, for 5-10 frequency bands. Run the arithmetic and you get days-to-weeks of GPU cluster time per survey, and that is already after using every computational trick available. This section catalogues the tricks.

1. Source encoding, the central trick

The wave equation is linear in the source: if shot $i$ produces data $d\_i$, the combined source $\\sum\_i c\_i s\_i$ produces data $\\sum\_i c\_i d\_i$. Pick random signs $c\_i \\in \\{-1, +1\\}$, build a super-source and a super-data, and run one wave simulation that collectively informs every shot's gradient:

The first term is $\\sum\_i g\_i$, the true FWI gradient we wanted. The second term is cross-talk: contributions from shot $i$'s forward wavefield correlated with shot $j$'s adjoint wavefield. Because $c\_i c\_j$ averages to zero over many independent random encodings, the cross-talk averages out over iterations, provided fresh random codes are drawn at every iteration; a fixed encoding leaves the same cross-talk in every gradient. One forward and one adjoint solve per iteration instead of $N\_{\\text{shots}}$ of each: about $N\_{\\text{shots}}/3$ times fewer solves once the roughly 3× more iterations the cross-talk noise costs are paid (Krebs et al. 2009 report a few-fold increase).

2. What encoding buys, and what it costs

Figure 6.3 prices one inversion three ways, counting two wave solves (a forward and an adjoint) per shot per iteration per frequency band, and shows on a small line of eight shots what the encoded gradient looks like. Grow the survey with $N$, change the iteration penalty encoding pays, average the cross-talk over iterations with fresh or frozen codes, and switch the receivers to a towed streamer.

At the starting setting, 1000 shots and 50 base iterations, the naive inversion needs 100 000 wave solves. Encoding needs 300: one super-shot, two solves, 150 iterations, 333 times fewer. A mini-batch of 100 random shots per iteration needs 14 600 (73 iterations under the figure's assumed noise penalty), 6.8 times fewer. The price of encoding is in (d): on the figure's line of eight shots, one encoded gradient carries cross-talk 180 % the size of the true gradient in (c), a level that grows with the number of shots blended into a super-shot. Averaged over 100 iterations with fresh codes it falls to 20 %, roughly as $1/\\sqrt{M}$ (f); with the codes frozen it never falls at all.

The saving grows with $N\_{\\text{shots}}$: at $N = 10\\,000$ it is $10\\,000 \\times 50 / 150 \\approx 3\\,333\\times$ in wave solves. The wall clock falls less, as the table in (h) shows: the encoded iterations run one after another, and a single super-shot cannot be spread over the shots of a cluster, so on 100 GPUs the mini-batch finishes first even though it needs 49 times more solves.

3. When encoded FWI breaks

• Moving-receiver (streamer) acquisition. Forming $d\_{\\text{enc}} = \\sum\_i c\_i d\_i$ requires every receiver to have recorded every shot. In towed-streamer data each shot has its own receivers, so the encoded observed data do not exist (in Figure 6.3 more than half of the super-shot's simulated data has no recorded counterpart); encoding is used mainly on fixed-spread land and ocean-bottom data, and streamer FWI relies on mini-batches instead.
• Strong amplitude variations between shots. A shot 10× stronger than the rest dominates the true gradient $\\sum\_i g\_i$ in conventional FWI too; what encoding adds is cross-talk. The strong shot takes part in every cross term $X\_{ij}$ it shares, so its cross-talk swamps the weak shots’ contributions, and the inversion needs more iterations, or amplitude balancing of the shots before encoding.
• Noise in the data. Encoding does nothing to suppress noise in the recorded data (swell, electrical hum): it enters the encoded residual as $\\sum\_i c\_i n\_i$ and must be removed in preprocessing or tolerated by a robust misfit, exactly as in conventional FWI.
• Salt imaging. At sharp reflectors, the cross-talk can add coherent artefacts that never fully wash out. Hybrid strategies (encoded outer, naive inner) are used.

4. Memory, storing the forward wavefield

The adjoint-state gradient requires the forward source wavefield $U\_s(x, z, t)$ at every grid point and every time step. For a 3D volume with $400 \\times 400 \\times 320 \\times 4000$ samples (single precision), that is ~800 GB per shot. Three standard mitigations:

• Checkpointing: save $U\_s$ at every $k$-th time step only. To compute the gradient at the skipped time steps, re-propagate forward from the nearest checkpoint. Memory drops by roughly $k$ at the cost of about one extra forward simulation, since each segment is recomputed once. When memory is fixed, Griewank's binomial (revolve) checkpointing gives the optimal schedule, with recomputation growing only logarithmically in the number of time steps.
• Random boundaries: instead of absorbing-boundary wavefield storage, randomise the velocity at the boundaries so outgoing waves return scrambled. Because nothing is absorbed, the wave equation stays time-reversible, so back-propagating from the final two snapshots reconstructs $U\_s$ exactly; the random boundary only scatters the boundary reflections so they do not correlate coherently in the gradient. Storage drops to two snapshots, at the cost of one extra propagation.
• Boundary-saving reconstruction: store the wavefield on a thin strip at the model boundary at every time step plus the final snapshot, then run the wave equation backward in time while re-injecting the saved boundary values. Storage is the strip rather than the volume, and the cost is one extra propagation.

5. Parallelism

• Shot-parallel: each shot's simulation is independent; distribute across GPU nodes. Embarrassingly parallel up to $N\_{\\text{shots}}$ workers; an encoded super-shot is one shot, so it gains nothing here.
• Domain-decomposed: split the model grid across GPUs within a node; each GPU handles its slab. Adds inter-GPU communication at slab boundaries but necessary for very large models.
• Frequency-parallel: in frequency-domain FWI the frequencies inside one group are independent and can be solved on separate GPU pools; across groups the multiscale strategy stays sequential, because each band starts from the previous band's model.

6. Hessian, or: why L-BFGS is the default

Pure steepest-descent FWI converges slowly because $J(m)$ has very different curvature along different model directions. Newton's method $m \\leftarrow m - H^{-1} g$ fixes this but requires the full Hessian $H = \\partial^2 J / \\partial m^2$, which needs $N\_{\\text{model}}$ Hessian-vector products (each about two extra wave solves per shot) to form and $N\_{\\text{model}}^2$ numbers to store, infeasible for $N\_{\\text{model}} \\sim 10^8$.

L-BFGS (limited-memory BFGS) approximates $H^{-1}$ from the last $m$ gradient-model pairs ($m = 5$ to $20$), giving near-Newton convergence at storage cost $2m\\,N\_{\\text{model}}$. Combined with an approximate diagonal pre-conditioner (Gauss-Newton on the diagonal, or an illumination-based pseudo-Hessian), L-BFGS is the most common default solver in production packages; nonlinear conjugate gradient and truncated Gauss-Newton are the main alternatives. Convergence is typically 3-5× faster than steepest descent per outer iteration, paying for itself immediately.

7. A realistic computational budget

For a 3D deep-water sub-salt project, 10 km × 15 km × 8 km at $\\Delta x = 25\\ \\text{m}$, 8 Hz, 5000 shots:

• Naive: 2 solves × 5000 shots × 200 iters × 6 bands = 12 million solves ⇒ about 2.3 years on a 100-GPU cluster at 10 GPU-minutes per solve.
• Encoded: 2 solves × 1 super-shot × 600 iters × 6 bands = 7200 solves ⇒ 1200 GPU-hours, about 12 hours of the same cluster's capacity. But the 3600 iterations run one after another and one super-shot cannot be split across shots, so without domain decomposition the wall clock is about 50 days on one GPU at a time.
• Mini-batch $k = 500$: 2 solves × 500 shots × 219 iters × 6 bands = 1.3 M solves ⇒ about 91 days, with 5 GPU batches per iteration.

Encoded FWI is the only option that fits a small budget in GPU-hours when the acquisition allows it (fixed-spread land or ocean-bottom nodes); splitting each super-shot over many GPUs, or running several super-shots per iteration, turns that saving into wall clock. For streamer data, mini-batches and shot-parallel hardware carry the load. Set Figure 6.3 to 5000 shots, 200 base iterations, 6 bands and a batch of 500 to reproduce these numbers.

Where this goes next

Section 6.4 moves from acoustic to elastic/anisotropic physics, what happens when the acoustic wave equation is wrong and converted waves or anisotropy-induced travel-time errors dominate the residual.

References

• Virieux, J., Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1.
• Pratt, R. G. (1999). Seismic waveform inversion in the frequency domain, Part 1. Geophysics, 64, 888.
• Etgen, J., Gray, S. H., Zhang, Y. (2009). An overview of depth imaging in exploration geophysics. Geophysics, 74, WCA5.
• Tarantola, A. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259.
• Krebs, J. R., Anderson, J. E., Hinkley, D., et al. (2009). Fast full-wavefield seismic inversion using encoded sources. Geophysics, 74, WCC177.