Source-linked AI summary
Application of Optimal Transport and the Quadratic Wasserstein Metric to Full-Waveform Inversion
Yunan Yang, Björn Engquist, Junzhe Sun, Brittany D. Froese
TL;DR
FWI’s standard L2 misfit is vulnerable to cycle skipping and noise, motivating an alternative that better handles waveform differences. The paper couples quadratic Wasserstein misfits to adjoint-state FWI, tests trace-by-trace and global formulations on three 2D models, and reports effective cycle-skipping mitigation with practical trace-by-trace cost. The method remains constrained by normalization-related nonconvexity and trace-by-trace scaling artifacts.
Problem
The commonly used L2 misfit in FWI suffers from cycle skipping and sensitivity to noise, limiting robust waveform-based inversion.
Method
The paper embeds the quadratic Wasserstein metric into adjoint-state FWI and compares trace-by-trace and global W2 formulations on three synthetic 2D models.
Results
W2-based FWI alleviates cycle skipping in numerical examples, while trace-by-trace W2 requires less than 1.1 times L2’s runtime and global W2 takes 3 to 4 times L2’s runtime.
Takeaways & Limitations
Quadratic Wasserstein misfits are a promising candidate for seismic inversion, with trace-by-trace W2 offering computation comparable to standard L2 FWI.
Takeaways & Limitations
Linear normalization is not convex with respect to simple shifts, and trace-by-trace rescaling can produce nonphysical adjoint-source variations and non-uniform data-misfit contributions.
Abstract
from arXiv · showhide
Conventional full-waveform inversion (FWI) using the least-squares norm ($L^2$) as a misfit function is known to suffer from cycle skipping. This increases the risk of computing a local rather than the global minimum of the misfit. In our previous work, we proposed the quadratic Wasserstein metric ($W_2$) as a new misfit function for FWI. The $W_2$ metric has been proved to have many ideal properties with regards to convexity and insensitivity to noise. When the observed and predicted seismic data are regarded as two density functions, the quadratic Wasserstein metric corresponds to the optimal cost of rearranging one density into the other, where the transportation cost is quadratic in distance. The difficulty of transforming seismic signals into nonnegative density functions is discussed. Unlike the $L^2$ norm, $W_2$ measures not only amplitude differences, but also global phase shifts, which helps to avoid cycle skipping issues. In this work, we build on our earlier method to cover more realistic high-resolution applications by embedding the $W_2$ technique into the framework of the adjoint-state method and applying it to seismic relevant 2D examples: the Camembert, the Marmousi, and the 2004 BP models. We propose a new way of using the $W_2$ metric trace-by-trace in FWI and compare it to global $W_2$ via the solution of the Monge-Ampère equation. With corresponding adjoint source, the velocity model can be updated using the l-BFGS method. Numerical results show the effectiveness of $W_2$ for alleviating cycle skipping issues and sensitivity to noise. Both mathematical theory and numerical examples demonstrate that the quadratic Wasserstein metric is a good candidate for a misfit function in seismic inversion.
1. Introduction
FWI seeks high-resolution subsurface models by iteratively matching synthetic and recorded seismic data, but conventional L2 misfits can suffer from cycle skipping and noise sensitivity. This paper develops W2-based alternatives grounded in optimal transport and applies them to realistic 2D FWI problems.
- FWI iteratively updates a subsurface model and synthetic data to reduce the misfit with recorded seismic data.
- The objective is to match synthetic and recorded waveforms comprehensively so that all waveform information contributes to the data misfit.
- The widely used L2 misfit suffers from cycle skipping and sensitivity to noise, while FWI also faces the inverse problem’s ill-posedness.
- Optimal transport compares seismic signals as density functions by finding the lowest-cost mass-transport plan, with W2 using quadratic transportation costs.
- Because seismic signals do not naturally satisfy nonnegativity and equal-mass constraints, they require normalization; earlier positive-negative separation is less effective for larger adjoint-state problems.
- The paper embeds W2 into adjoint-state FWI, introduces trace-by-trace and global comparisons, and evaluates them on Camembert, Marmousi, and 2004 BP models.
2. Theory
This section formulates FWI misfits and develops quadratic Wasserstein optimal transport for seismic data, including its one-dimensional and higher-dimensional constructions. It highlights W2's convexity for relevant model variations and reduced sensitivity to noise.
- FWI formulation: W2 can be applied trace-by-trace in one dimension or globally to the full data set, producing different misfit and adjoint-source formulations.The trace-by-trace formulation sums one-dimensional comparisons across the total number of traces, while the global formulation compares the complete data sets.
- FWI formulation: FWI minimizes a misfit between simulated and observed data to recover unknown subsurface model parameters.The model is updated through a gradient-based iterative scheme for a PDE-constrained optimization problem.
- Optimal transport formulation: The quadratic Wasserstein metric compares seismic signals as equal-mass nonnegative densities by minimizing quadratic-cost transport between them.The optimal map rearranges one distribution into the other, and its minimum transport cost defines the metric.
- One-dimensional transport: For one-dimensional densities, the optimal map is the unique monotone rearrangement T = G^-1 ◦ F, computed from cumulative distribution functions and their inverses.Signals must be normalized so they are positive, supported on [0, 1], and have total mass 1.
- Higher-dimensional transport: In higher dimensions, the optimal map is obtained through cyclical monotonicity and the Monge-Ampère equation with a convex potential.The substitution T(x) = ∇u(x) converts the transport constraint into the Monge-Ampère equation, after which the squared metric is computed from the map.
- Convexity and noise: The squared W2 metric is convex for data shifts, dilation, and partial amplitude changes, and is substantially less sensitive to noise than the traditional L2 norm.For typical seismic data, the effect of noise is expected to be negligible, including when noise has order-one amplitude with many data points.
3. Numerical scheme
The numerical scheme preprocesses seismic signals for Wasserstein comparison, computes trace-wise or global W2 misfits, and derives gradients for adjoint-state inversion.
- 3.1. Data normalization.: Signals are shifted and rescaled into positive, equal-mass data before applying Wasserstein-based FWI.The linear transformation selects c so f + c and g + c are positive, then rescales both signals to a common total mass.
- 3.1. Data normalization.: Linear normalization preserves waveform phase and supports differentiation, but does not preserve convexity with respect to simple shifts.The normalization maintains local extrema and has regularity favorable to the adjoint-state method, while its shift convexity remains unresolved.
- 3.2. Compare trace by trace: W2: The trace-by-trace method computes one-dimensional W2 values from cumulative distributions and their inverses, then sums them over receivers.Discrete cumulative distributions are obtained by numerical integration, with inverse evaluation accelerated by binary search and interpolation.
- 3.2.2. Computation of adjoint source.: The trace-wise misfit has O(N log(N)) complexity and its Fréchet gradient supplies the adjoint source for inversion.The gradient is derived from the first variation of the squared Wasserstein metric using finite-difference operators and the lower triangular integration matrix.
- 3.3. Compare global data: W2: Global W2 compares the entire data set through a higher-dimensional Monge–Ampère equation rather than a one-dimensional exact formula.The equation is discretized with finite differences, including an almost-monotone formulation that enforces convexity and is solved using Newton’s method.
- 3.3.2. Computation of adjoint source.: For global W2, linearizing the discrete Monge–Ampère equation yields a gradient formulation that reuses its formal Jacobian.The Jacobian is already inverted during Newton iterations, linking the Monge–Ampère solve to adjoint-source computation.
4. Computational results
The experiments compare L2, trace-by-trace W2, and global W2 on synthetic seismic models. W2-based methods address phase-related cycle skipping, with global W2 producing more accurate Camembert results than trace-by-trace W2.
- Experimental design: The computational study applies W2 both trace by trace and to entire data sets, comparing both approaches with conventional L2 inversion.Global W2 is restricted to smaller-scale models because of limitations in the current Monge-Ampère solver.
- 1D case study: W2 adjoint sources shift mass to correct phase differences, whereas L2 primarily corrects amplitude differences associated with cycle skipping.The W2 adjoint source is smoother and has no DC component, supporting numerical stability in quasi-Newton updates.
- Camembert model: Both trace-by-trace and global W2 converge in 10 l-BFGS iterations on the Camembert model.The Camembert inversion starts from a homogeneous 3km/s model with a 3.6km/s circular inclusion.
- Camembert model: L2 converges to a local minimum after 100 iterations, while W2 adjoint sources provide an update direction that avoids the reported cycle-skipping behavior.The L2 adjoint source has positive-negative-positive components, unlike the negative-positive W2 pattern.
- Camembert model: Global W2 produces a more accurate Camembert velocity model than trace-by-trace W2 in terms of numerical velocity error.The global and trace-by-trace approaches use different formulations of the misfit and adjoint source.
4.3. Scaled Marmousi model with global W2 misfit computation.
On a scaled Marmousi model, global W2 avoids the local-minimum behavior observed with conventional L2 from a highly smoothed initial model.
- Scaled Marmousi setup: Global W2 avoids the local minimum encountered by conventional L2 in the scaled Marmousi inversion.The initial model is a heavily Gaussian-smoothed version of the true velocity model.
- Scaled Marmousi comparison: L2 produces spurious high-frequency artifacts attributed to its point-by-point amplitude comparison.The comparison uses global W2 and conventional L2 misfit functions.
4.4. True Marmousi model with trace-by-trace W2 misfit computation.
Trace-by-trace W2 recovers Marmousi structure more effectively than L2 from a highly smoothed initial model and reaches a relative misfit of 0.1 in 20 iterations.
- Experimental setup: The trace-by-trace experiment uses 5Hz Ricker sources, 30m spatial and 30ms temporal discretization, and removes frequencies from 0 to 2Hz.The true Marmousi model is 3km deep and 9km wide.
- Trace-by-trace computation: Trace-by-trace W2 normalizes each receiver trace before solving the one-dimensional optimal-transport problem and combines the resulting derivatives into an adjoint source.The adjoint source is propagated backward to generate the velocity gradient.
- Gradient comparison: In the first iteration, W2 focuses its gradient on the Marmousi model’s peak, while the L2 gradient remains comparatively uniform.The W2 gradient’s darker areas match multiple features of the true velocity model.
- Inversion results: After 300 l-BFGS iterations, trace-by-trace W2 correctly inverts most Marmousi details, whereas L2 retains spurious high-frequency artifacts.The inversion starts from a true model smoothed with a Gaussian filter of deviation 40.
- Convergence: W2 reduces the relative misfit to 0.1 in 20 iterations, while L2 converges slowly to a local minimum.The convergence behavior is shown in the comparison curves for the two inversion methods.
4.5. Inversion with the noisy data.
Trace-by-trace W2 remains effective with strongly noisy Marmousi data and outperforms a later L2 refinement from the same improved starting model.
- BP model: The BP experiment also compares L2 and global W2 inversion results using true and initial velocity models.The supplied figure passages identify the BP true and initial models and the corresponding inversion results.
- Noisy-data inversion: Even when noise power exceeds signal power, W2 converges reasonably well and still recovers most Marmousi features.The noisy-data experiment uses uniform random iid noise with an SNR of −3.4716 dB.
- W2 initialization versus L2: The W2-initialized L2 result has 0.4211 relative error, compared with 0.2664 for trace-by-trace W2 after comparable subsequent iterations.Both results recover most Marmousi features, but the L2 approach converges more slowly even from a good initial model.
4.7. 2004 BP Model with global W2 misfit computation.
Global W2 inversion on the modified 2004 BP model recovered the top salt more reasonably than L2, which converged to a cycle-skipped low-velocity anomaly beneath it.
- The experiment compares global W2 and conventional L2 on a modified BP 2004 model with complex deep-water Gulf of Mexico geology.The initial model contains a smoothed background without the salt.
- Global W2 computes its misfit by numerically solving the Monge-Ampère equation.The inversions were stopped after 100 l-BFGS iterations.
- W2 recovered the top salt reasonably well, whereas L2 converged to a low-velocity anomaly immediately beneath the top salt.The anomaly is identified as typical of cycle skipping in FWI.
4.8. 2004 BP Model with trace-by-trace W2 misfit computation.
Trace-by-trace W2 produced informative early gradients, recovered BP salt-body shapes, and reduced the relative misfit rapidly compared with L2.
- The figures report the final data residual and first-iteration gradients for trace-by-trace W2 and L2.
- In the first iteration, W2 concentrated the inversion on the upper salt, while the L2 gradient was not very informative.The W2 gradient’s darker area matched the salt region in the velocity model.
- After 300 l-BFGS iterations, trace-by-trace W2 constructed the salt-body shapes, whereas L2 failed to recover their boundaries.
- Trace-by-trace W2 reduced the relative misfit to 0.1 in 20 iterations, while L2 converged slowly to a local minimum.The convergence comparison is shown in Figure 26.
5. Discussion on two ways of using W2
Trace-by-trace W2 is substantially cheaper than global W2 and performed effectively in the reported BP and Marmousi experiments, but scaling and numerical artifacts remain important considerations.
- Trace-by-trace W2 requires less than 1.1 times the runtime of L2, whereas global W2 took 3 to 4 times the L2 runtime.The cost difference reflects one-dimensional optimal transport versus solving a two-dimensional Monge-Ampère problem.
- Trace-by-trace rescaling can create nonphysical amplitude variations and non-uniform backgrounds in the adjoint source.These effects may produce non-uniform contributions of data misfits to velocity updates, motivating more careful scaling.
- Global Monge-Ampère computation requires enough data points to resolve steep gradients because its discretization error depends on the target profile’s Lipschitz constant.Insufficient resolution effectively regularizes the data before solving the equation.
- Trace-by-trace optimal transport uses exact one-dimensional formulas and remains accurate for highly non-smooth data.
- W2’s noise insensitivity can coexist with oscillatory artifacts, which may arise from the numerical PDE solution and the metric’s noise insensitivity.The global-W2 Camembert result had a model error 25% smaller than the trace-by-trace result.
- Additional L2 iterations did not obviously improve the trace-by-trace W2 velocity resolution, while 20 W2 iterations produced a lower model error than L2 after the same iteration count.The W2 result can also serve as an initial model for higher-frequency L2 FWI.
6. Conclusion
The paper develops high-resolution FWI using W2 with adjoint-state optimization and demonstrates successful applications across several 2D models. It reports comparable speed to L2 with greater accuracy, faster convergence, and reduced cycle skipping, while identifying normalization as an area for improvement.
- The proposed technique couples the W2 misfit with efficient adjoint-source computation for high-resolution FWI.It is applied successfully to the Marmousi, 2004 BP, and Camembert models.
- Global W2 uses a Monge-Ampère equation, while trace-by-trace W2 provides comparable inversion results.
- W2-based FWI is reported as as fast as standard L2 FWI, but more accurate, faster-converging, and less affected by cycle skipping.
- Signal scaling and normalization strongly affect W2 results, and the linear normalization that worked best for large-scale inversion does not satisfy the theoretical shift-convexity requirement.The paper identifies improved normalization as a future research direction.
Appendix A. Derivation of Equation (33)
The appendix derives the Fréchet derivative of the quadratic Wasserstein expression with respect to a continuous density function, then discretizes the result for the numerical scheme.
- The derivation assumes continuous density functions f(t) and g(t) on the interval [0, T0].
- The analysis perturbs f by δf and examines the resulting variation of the Wasserstein expression as a functional of f.
- A Taylor expansion of the monotone inverse function G−1 is used to obtain the first variation.
- The first variation yields the Fréchet derivative of the Wasserstein expression with respect to f.
- The numerical scheme discretizes this derivative and derives equations (31) and (33).