In classical waveform inversion (Lailly, 1983; Tarantola, 1984; Virieux and Operto, 2009), we estimate the earth parameter, m, by minimizing the data residuals. In a time domain formulation, the classic objective function reads J m = 1 2 ∑ s , r ∫ 0 T W s , r t p s , r t − d s , r t 2 d t , (1) where s is the source index, r is the receiver index, p _s,r is the synthetic trace, d _s,r is the observed data trace, and W _s,r is a time window to potentially select the events of interest. The recording length is T. With long-offset data, this window can, for instance, be used to select only the early arrivals, as explained later. The adjoint source is the partial derivative of the right-hand side in Eq. 1 with respect to the synthetic data: s t = W s , r 2 t p s , r t − d s , r t , (2) with s(t) as the adjoint source for the trace indexed by s, r (we omitted the lower index).

Our 2D acoustic forward wave simulator is based on a staggered-grid finite-difference scheme (Virieux, 1986), with a free surface boundary condition on top. FWI is based on a local optimization technique or gradient method to iteratively update the unknown parameters (model and source wavelets). The gradient is computed via the adjoint state method (Liu and Tromp, 2006; Plessix, 2006). Briefly speaking, it usually involves a couple of steps. First, a forward wave propagation is conducted for the given model. The border values and the synthetics are saved either in memory or a disk. Then, the synthetic data are compared to the observed data using a given objective function, producing the misfit and the adjoint source. Finally, the adjoint wave fields are simulated with the adjoint source using an adjoint simulator, while, at the same time, the forward wavefields are reconstructed from the saved border values. The accumulated over time product of the adjoint and the reconstructed forward wavefields provides the gradient.

The gradients, which usually contain numerical artifacts for many reasons, are smoothed with an elliptic Gaussian filter. At the first scale, the standard deviation of the Gaussian smoothing filter is 750 m and the vertical one is 150 m. At the final scale, the standard deviation is 150 m and the vertical one is 30 m. To compensate for the decay of the gradient energy with depth, we multiply the gradient by z + 300 , where z is the depth in meter. Finally, the minimization of the objective function is performed using an l-BFGS quasi-Newton algorithm (Nocedal, 1980).

For the seismic trace indexed by s, r, we denote p ̄ ( t ) = W s , r ( t ) p s , r ( t ) and d ̄ ( t ) = W s , r ( t ) d s , r ( t ) . The cross-correlation function between the windowed seismic traces is defined as c τ = xcorr p ̄ , d ̄ = ∫ 0 T p ̄ t d ̄ t − τ d t . (3) The normalized cross-correlation function is defined based on the aforementioned cross correlation as J c = − ∫ − T T c 2 τ P τ d τ ∫ − T T c 2 τ d τ , (4) where P ( τ ) = exp ( − τ 2 2 ζ 2 ) is a Gaussian penalty function. ζ is a user-defined parameter that controls the width (proportional to the standard deviation) of the penalty function. We choose ζ to be 5% of the whole recording length ζ = 0.05T in the following inversions. The minus sign is to transform the equation into a minimization problem. The adjoint source is s c t = W s , r t ∫ − T T c 2 τ d τ ∫ − T T P τ − 2 J c c τ d ̄ t − τ d τ . (5)

Because of the adjoint state method (Liu and Tromp, 2006; Plessix, 2006), it is very convenient to switch from L2 FWI to FWI with the normalized cross-correlation J _c . In an existing FWI program, only the misfit value and the adjoint source need to be changed.

Both L2 and the normalized cross-correlation function J _c are 1D comparisons of the seismic data along the trace. L2 is considered a pixel-by-pixel comparison. The optimal transport, depending on the type, could be 1D or 2D. Its 2D version could exploit the lateral coherency along the receiver axis besides the coherency along the time axis. In this case, the data difference for a shot indexed with s is J t = max φ s ∈ B L i p 1 ∑ r ∫ 0 T φ s x r , t W s , r t δ p s , r t d t Δ x r , (6) where Δx _r is the receiver spacing, W _s,r (t) is a weighting function introduced to represent the time-windowing technique, and the signed residual δp _s,r = p _s,r − d _s,r is the difference between synthetic and observed shot gathers. BLip₁ is the space of bounded 1-Lipschitz functions defined in the (x _r , t) shot-gather space (Métivier et al., 2016a; Poncet et al., 2018) by ∀ x r , t : | φ s x r , t | < λ , ∀ x r , t : | φ s x r + δ x r , t − φ s x r , t | < δ x r , | φ s x r , t + δ t − φ s x r , t | < δ t , (7) where λ is a pre-defined constant, δx _r is any positive increment in the spatial direction, and δt is any positive temporal increment. The first inequality ensures that the amplitudes of φ _s are bounded; the last two inequalities prevent abrupt variations. The choice of λ was discussed by He et al. (2018). For the 1D optimal transport function, the constraint along the spatial axis is removed. The maximization problem (Eq. 6) can be solved efficiently through proximal splitting techniques (Métivier et al., 2016b). Once the optimal φ ̄ s ( x r , t ) is determined, the corresponding adjoint source is simply s t t = φ ̄ s x r , t W s , r t . (8)

This is a 2D algorithm, and we repeat it for each shot panel in our following inversion examples. Again, thanks to the adjoint state method (Liu and Tromp, 2006; Plessix, 2006), switching from a least-squares objective function to the OT objective function in an existing FWI code only requires changes in the misfit value and, specifically, in the adjoint source.

A quick way to evaluate an objective function is to observe its behavior with moving Ricker wavelets. Figure 1A shows such a test. A 3-Hz Ricker wavelet centered in the middle is taken as an observed seismic trace. The synthetic trace is a shifted 3-Hz Ricker wavelet. The L2 norm objective function shows two local minima; the valley near the true position is really sharp. This is consistent with the fact that cycle skipping is an issue for the L2 norm, but it can help FWI converge fast and produce a high-resolution model if we start with a model within the basin of attraction. For the normalized cross-correlation function J _c , we choose the penalty function P(τ) with ζ = 0.3T = 1.2 s. In Figure 1A, we re-scaled the range to [0,1]. The misfit curve is similar to the travel time-based objective function. There are no local minima. However, the wide valley suggests that FWI, in this case, will converge to a smoother model. We also show the misfit curve with optimal transport function J _t with the 1D algorithm. There are no local minima, and the basin of attraction is a little bit wider than that of the L2 function. However, if the synthetics is too far from the observed trace, the misfit J _t is flat, indicating that there is no sensitivity to the variation in synthetics (Figure 1B). The adjoint sources under J _c and J _t are also shown (the L2 adjoint source is the direct difference between the synthetics and observed trace; thus, it is omitted). For the J _c function, the position of the adjoint source signal aligns with the synthetic trace. This behavior is similar to that of the travel time-based objective function (Dahlen et al., 2000; Nolet, 2008), where the adjoint source is proportional to the derivative of the synthetics p ̇ s , r ( t ) . For the J _t function, the imprints of the synthetic and observed trace are kept in the adjoint source. However, the adjoint source is less oscillating than the synthetic and the observed trace. This is consistent with the fact that the misfit curve of J _t is more convex than that of the L2 function (Figure 3A).

We shift the phase of the synthetics by 90° and calculate the misfit curves and the adjoint sources for these three objective functions in Figure 2. For both the L2 function and J _t function, there are local minima. At the correct travel time t = 0, the misfit values are neither local minima nor local maxima. However, for the J _c objective function, the misfit curve is the same as that shown in Figure 2B. Both of the adjoint sources of J _c and J _t change. The phase shift in the synthetic trace (Figure 2A) is taken into account by the J _c function.

A realistic model mimicking the geology of the Middle East area is taken as the true model. The model contains a succession of shale–sand layers with velocities between 2 and 3.5 km/s, depending on the depth, and carbonate layers with velocities between 3.5 and 6 km/s (Figure 3). Most of the layers have a thickness of several hundred meters.

The acquisition mimics a standard fixed-spread vibroseis acquisition. The observed dataset contains 48 sources, distributed every 300 m, and 300 receivers, spaced every 50 m on the free surface. The observed data are modeled using an acoustic forward operator with a free surface boundary condition. The source wavelet is a 3-Hz Ricker wavelet. The total recording time is 8.192 s. A 1D linearly increasing velocity from 2.2 km/s at the top to 5.5 km/s at the bottom is taken as the initial v _p model (Figure 4). The initial density is linked to the initial v _p model via Gardner’s law (Gardner et al., 1974; Hamilton, 1978). The initial density is overestimated in the deep area.

The observed data are shown in Figure 5. The complex observed seismogram contains various waves. Diving waves penetrate deep. Near-offset reflections and inter-bed reflections are due to the lateral layers in the top area. Scattered events are due to the lateral curvature in the middle of the model shown in Figure 3.

The conventional inversion is to start from low frequency. We consider the L2 function. The first frequency band is [0.6, 1.5] Hz, which is sufficiently low in typical exploration geophysics. The density is re-calculated from v _p at each iteration via an empirical relationship. By passively updating the density parameter, structural information is introduced from the inverted v _p model to density.

The linearly increasing initial model is too far from the true model. Therefore, when estimating the source wavelets in the initial model, the initial estimated source wavelet is completely different from the true source wavelet: its amplitude is 1,000 times smaller, although the frequency band is more or less consistent. Consequently, l-BFGS only updated the model for two iterations before having difficulty determining a suitable step length. The final model is almost identical to the initial model. We do not move to a higher-frequency band. Our rationale is that if the inversion is trapped in a local minimum at low frequency, moving to high frequency induces more local minima.

We replace the L2 function with OT objective function and repeat the inversion in this low-frequency band [0.6, 1.5] Hz. Nevertheless, a similar phenomenon is observed, probably because the estimated source wavelet is very far from the true source wavelet.

Instead of a frequency continuation strategy, we can also consider an alternative inversion strategy based on seismic events. We do not need accurate windows to isolate seismic events. Our idea is to match the early arrivals first. When the model is better recovered, we gradually increase the time windows to include more and more complex events. The requirement on the design of time windows is that it should be smooth, such that later events are smoothly included. We consider Gaussian time windows. For each trace, a Gaussian window is W s , r t = exp − t − t s , r 0 2 2 σ 2 , (9) where t ₀ is around the first break of each trace. Here, the frequency continuation inversion strategy is not used. The inversion strategy depends on the standard deviation σ. For the total recording length T = 8.192 s, the first three inversion stages are conducted with the ratios σ T = 0.2 , 0.3 , 0.5 , while in the fourth stage, the entire seismograms are included without using time windows. During the inversion, density is passively updated.

The inversion is much better than the one with the frequency continuation strategy. Figure 6 shows the final recovered source wavelets. Both L2 and OT inversions recover wavelets close to the true source wavelet. Figure 7 shows the final recovered v _p velocity models. Some structures are recovered, but there are artifacts in both inversions.

We observed that an event-based inversion strategy helps recover the structure, but it is not sufficient with L2 or the OT objective function. We propose to apply the normalized cross-correlation function J _c first and then conduct an OT inversion. We still consider the event-based inversion strategy. In the J _c inversion, we select a long offset trace as the source wavelet and keep it fixed during the inversion. When the OT inversion is conducted, both the source wavelets and the model parameters are simultaneously updated.

The normalized cross-correlation function J _c is less sensitive to source wavelets, as shown in Figure 1. We choose a long offset seismic trace at X=15 km (blue line in Figure 5). We calculate the envelope of the whole trace and cut the seismic trace defined by the first envelope wavelet (we select the seismic trace, not the envelope). Figure 8 shows how we select the seismic wavelet. Its Fourier transform shows that the selected seismic wavelet has a different frequency content compared to the true source wavelet. Figure 9 shows the synthetic data corresponding to the cut seismic trace in the linearly increasing initial model. The initial seismograms are quite simple.

The cross-correlation could map multiple events as sidelobes. This affects the adjoint source and the consequent gradient, leading to artifacts in the model. It is important to not include too many events under J _c . Otherwise, it is difficult to obtain an acceptable inversion. Figure 10 shows an inversion under J _c without using any time window. The initial source wavelet is the cut trace shown in Figure 8. The gradient does not show a correct direction update between the first and third dashed lines. The final model is not close to the true model.

The inversion using J _c is much better when the Gaussian time windows are used than the one without windows. For the total recording length T = 8.192 s, four inversion stages are conducted with the ratios σ T = 0.02 , 0.05 , 0.1 , 0.2 . The first Gaussian time window keeps only the early arrival events. The gradient is shown in Figure 11. In comparison, we show the gradient of WET where a true source wavelet is used. The similarity indicates that J _c combined with Gaussian time windows can produce correct gradients even with a wrong source wavelet. The wider windows allow us to include more complex events. We stop at the ratio of 0.2. Later stages with bigger ratios introduce too much cross-talk between events to the adjoint source, and the model updating becomes negligible. A shot gather from the final synthetic data shown in Figure 12 shows that the global pattern of the observed data is captured. The final velocity model v _p is shown in Figure 13. The long-wavelength structures are retrieved. The fine-scale structures are not very satisfactory because of the cross-talk between seismic events when the Gaussian window becomes wide.

We continue the inversion by following an OT objective function. The Gaussian time windows are defined with ratios σ T = 0.2 , 0.3 , 0.5 in three stages. Figure 14 shows the estimated source wavelets for these three stages. The first wavelet clearly deviates from the true source wavelet. However, later source wavelets become closer and closer to the true source wavelet. The final v _p and ρ are shown in Figure 15. The continuous thin-layered structures along the first dashed line are well-retrieved in the v _p model. The inverted v _p model is much cleaner than the starting model and shows high-resolution, high-fidelity structures. The inverted ρ is overestimated in the deeper area below the second dashed line. Above the second dashed line, thin-layered structures are obtained. A shot gather from the final synthetic data is shown in Figure 16A. The final residuals shown in Figure 16B indicate that the early arrivals are not perfectly matched. The trace-by-trace comparison shown in Figure 17 indicates that the phases are fitted, but the amplitudes are not perfectly matched. The reason is that the Gardner law overestimated the density in the deep area; therefore, the inversion does not perfectly fit the amplitudes of the early arrivals.

We observed, in the Middle East model, the presence of abundant diving waves in the data. To show the generalization of our proposed inversion flow, we now use the SEAM model, which is shown in Figure 18. Geologic layers are generally not horizontal for the SEAM model but are folding. There is also displacement of geologic layers. The velocity v _p increases slowly from top to bottom, unlike in the Middle East model. The initial velocity model v _p is also a 1D linearly increasing model, and the density is calculated from v _p via the same empirical relationship as that used in the Middle East model inversion. Figure 19 shows the initial v _p and ρ. We also overlap the traced rays for one shot in the initial model. The deepest ray penetrates only to approximately 1 km. The difference between the SEAM model and the Middle East model can also be observed in the seismograms. The observed data shown in Figure 21A are dominated by reflections.

The inversion setting is the same as in the Middle East model inversion. We use the same source wavelet, consider a fixed seismic acquisition, and also consider the same inversion strategy. As we did in the Middle East model inversion, we select a long offset seismic trace and cut the first wavelet as the initial source wavelet. The final model after J _c inversion is shown in Figures 20A, B, which shows that the model is largely updated. Then, we continue the inversion with the OT objective function J _t . The final model is shown in Figures 20C, D. The twisted geologic structures are well-recovered from top to bottom. The final synthetic is shown in Figure 21B, which captures the pattern of the observed data shown in Figure 21A. The residual shown in Figure 21C also confirms that most of the events are matched. In the true model, the structures in the leftmost part would reflect seismic waves outside of the acquisition. This explains that the model is not well-recovered in the leftmost part. This is also consistent with the fact that the non-zero part of the residual is mainly short offset reflections.

We apply our inversion workflow to a real-land seismic dataset recorded in Southwest China. The land seismic data have been preprocessed to remove the surface waves and to correct for varying topography. The seismic survey is along a straight line. It is 75-km long. The shot spacing is not regular: it varies from 20 m to 200 m. Totally, we have 619 shots. For each shot, the offset is from approximately −6 to 6 km, with an inter-receiver spacing of 20 m. Thus, each shot is recorded by less than 600 geophones. The recording time is 2 s, and the time interval is 1 m. The processed observed data are shown in Figure 22A, which contains direct waves and reflections. The maximum frequency is approximately 15 Hz, while the lowest frequency is approximately 5 Hz. The initial model, shown in Figure 23A, is provided by an industrial routine process including ray tracing tomography. The initial synthetics are shown in Figure 22B. The interleave comparison shown in Figure 22C shows the data fitting of initial synthetics to the observed data.

Following our inversion workflow, we apply the normalized cross-correlation objective function in the first step. The model is updated a little bit. Then, in the second step, we parallel-run two inversions with the conventional least-squares objective function and OT objective function. In this second step, we divide the frequency band into three stages with the maximum frequency at 7 Hz, 9 Hz, and 15 Hz. In each frequency band, we run two inversions with a Gaussian window σ = 0.2T and without windows. The conventional inversion wildly updates the model, as shown in Figure 23B. The least-squares function is a point-by-point comparison process. Thus, the pattern of seismic events is not honored. When the adjoint source is back-propagated to assemble the gradient, artifacts emerge. At the end, model updating becomes inaccurate. Figure 23C shows the final model recovery with the OT objective function. Coherent structures are revealed. Fine structures are reconstructed. A well log comparison at X = 45.5 km is shown in Figure 24. The velocity down to 3 km is overestimated because we use a constant density inversion, so the effects of density are not taken into account in the final v _p model. High-frequency model variation is recovered. The final synthetics, shown in Figure 25A, are closer to the observed data than the initial synthetics. The inter-leave comparison shown in Figure 25B shows the data fitting of final synthetics to the observed data.