# Geophysics

## ABSTRACT

In reflection seismology, full-waveform inversion (FWI) can generate high-wavenumber subsurface velocity models but often suffers from an objective function with local minima caused mainly by the absence of low frequencies in seismograms. These local minima cause cycle skipping when the low-wavenumber component in the initial velocity model for FWI is far from the true model. To avoid cycle skipping, we discovered a new wave-equation reflection traveltime inversion (WERTI) to update the low-wavenumber component of the velocity model, while using FWI to only update high-wavenumber details of the model. We implemented the low- and high-wavenumber inversions in an alternating way. In WERTI, we used dynamic image warping (DIW) to estimate the time shifts between recorded data and synthetic data. When compared with correlation-based techniques often used in traveltime estimation, DIW can avoid cycle skipping and estimate the time shifts accurately, even when shifts vary rapidly. Hence, by minimizing traveltime shifts estimated by dynamic warping, WERTI reduces errors in reflection traveltime inversion. Then, conventional FWI uses the low-wavenumber component estimated by WERTI as a new initial model and thereby refines the model with high-wavenumber details. The alternating combination of WERTI and FWI mitigates the velocity-depth ambiguity and can recover subsurface velocities using only high-frequency reflection data.

## INTRODUCTION

Traveltime inversion (Bishop et al., 1985; Stork, 1988; Luo and Schuster, 1991) and full-waveform inversion (FWI) (Tarantola, 1984; Pratt, 1999) are two different methods that geophysicists use to estimate subsurface velocity models from reflection seismic data. Both methods have advantages and disadvantages in estimating a subsurface model

Traveltime inversion minimizes only traveltime differences between observed (recorded) reflection data

Traveltime inversion has succeeded in applications to crosswell data (Luo and Schuster, 1991; van Leeuwen and Mulder, 2010) and refraction data (Zelt and Barton, 1998), where commonly only first arrivals are used to measure traveltime shifts. However, reflection traveltime inversion (RTI) is more difficult. One difficulty is that traveltime shifts vary with time. Manually picking arrivals in seismograms is the most reliable way to estimate these shifts, but it is time consuming. More automatic approaches compute crosscorrelations of seismograms to estimate time shifts (Luo and Schuster, 1991; Hale, 2009; van Leeuwen and Mulder, 2010). Unfortunately, when time shifts change rapidly (Hale, 2013), correlation-based methods still suffer from the cycle-skipping problem. A second difficulty in RTI is the ambiguity between the low-wavenumber background velocities and the depths of high-wavenumber reflectors (Stork and Clayton, 1986; Bube et al., 2002). This ambiguity can be resolved by the migration/demigration approach in tomography techniques (Lambaré, 2008; Woodward et al., 2008), which minimizes residual moveout (RMO) in migrated gathers. RMO-based tomography is not the focus of this paper, and we concentrate on only traveltime inversion and FWI, two wave-equation-based methods with objective functions defined in the domain of unmigrated data.

Conventional FWI does not differentiate

Difficulties in conventional reflection FWI are mainly due to local minima, which often arise when the initial background model is far away from the true model and the recorded data do not contain sufficient low frequencies. As a result, the recovery of high wavenumbers suffers from cycle skipping (Virieux and Operto, 2009), one primary cause of local minima. The other important reason for local minima is the nonlinear relationship between data and the model (Snieder, 1998). In the presence of cycle skipping, the high-wavenumber details

We propose a new wave-equation reflection traveltime inversion (WERTI) to update the low-wavenumber

In the following sections, we first introduce DIW to estimate traveltime shifts. Then, we analyze and compute the gradient in WERTI using the adjoint-state method (Plessix, 2006). We then test our method for a high-velocity-anomaly model and a part of the Sigsbee2A model (Paffenholz et al., 2002). Results demonstrate that by combining WERTI and FWI, we can use only high-frequency reflection data to recover low- and high-wavenumber components of velocity models.

## WERTI

Common methods for estimating time shifts are based on crosscorrelation (Luo and Schuster, 1991; van Leeuwen and Mulder, 2010). Crosscorrelation-based methods, however, may suffer from cycle skipping when shifts change rapidly. In this paper, we choose DIW to overcome the cycle-skipping problem in estimating time shifts. After posing an inverse problem for WERTI, we investigate the gradient for WERTI to update the low-wavenumber component of velocity models and propose an alternating inversion strategy.

### Traveltime shifts from dynamic warping

DIW (Hale, 2013) is an efficient extension to a classic dynamic warping algorithm for speech recognition (Sakoe and Chiba, 1978), and it can be considered a nonrigid image registration (Glasbey and Mardia, 1998; Zitová and Flusser, 2003). DIW obtains the traveltime shifts

Constraints in equation 6 control the rates at which shifts vary in directions

Figure 1 shows an example of dynamic warping applied to estimate shifts between two synthetic seismograms,

### Inverse problem

After we obtain the traveltime misfit

Various methods, including conjugate-gradient methods and quasi-Newton methods, can be used to solve the inverse problem implied by equation 8. All of the methods require computing the gradient of the objective function with respect to the model parameters.

### Gradient calculation

We expand equation 8 about some starting models

In this paper, we instead use the adjoint-state method (Plessix, 2006) to derive the gradients

The first term on the right side of equation 10 represents the zero-lag crosscorrelation of a forward source-side wavefield

With these virtual sources, the calculation of

### Comparison to FWI

FWI uses a different objective function that includes (directly) amplitude differences and (indirectly) traveltime differences:

As in the derivation above of gradients in WERTI, we can also get two gradients for FWI, with respect to

Meanwhile, in the context of conventional FWI, we aim to minimize the same data-difference objective function (equation 12). However, we do not distinguish gradients

### Which gradient?

The previous comparison highlights that WERTI and FWI can provide gradients with respect to the low-wavenumber background model

On the one hand, WERTI measures the kinematic difference that is more sensitive to

Although, as implied by equation 13, one can obtain for FWI a low-wavenumber gradient

## INVERSION PROCEDURE

Our inversion methodology is the following two-step combination of WERTI and FWI, both of which use a conjugate-gradient method:

begin with a smooth background, |

1) for do |

FWI (1 iteration): |

with |

obtain |

WERTI (1 iteration): |

with |

update |

update |

end for |

use |

2) for do |

iterate conventional FWI: |

compute gradient |

update |

end for |

In practice, inversion often begins with a smooth model

In the next iteration of step 1, we apply again one iteration of FWI after WERTI to obtain a new rough-model perturbation

True amplitude migration or least-squares migration can be used in step 1 to replace FWI and to approximate

When performing WERTI in step 1, we apply a spatial low-pass filter to the gradient

## EXAMPLES

To demonstrate the application of WERTI, we test the algorithm using a simple Gaussian high-velocity anomaly model and a part of the Sigsbee2A model. In both tests, we use only reflection data; we mute any direct and refracted arrivals.

### A simple velocity anomaly

We refer to the model

Figure 5a and 5b shows the FWI gradient

WERTI first synthesizes reflection data using the model shown in Figure 5b, and one of the common-shot gathers is displayed in Figure 6b. Figure 6c illustrates the difference between the true data (Figure 6a) and the synthetic data (Figure 6b) and highlights the cycle skipping, especially for offsets

More compelling evidence of cycle skipping is in the time shifts

After obtaining the time-varying time shifts

We can modify step 1 by replacing WERTI with FWI (equation 12) to minimize the data-difference objective function and then use the new low-wavenumber gradient

### Sigsbee model

Figure 10a displays the true velocity model that is part of the Sigsbee2A model created by SMAART JV (Paffenholz et al., 2002); Figure 10b shows an initial velocity, which increases as a function of depth. Except for the water layer, which is the same as in the true model, the initial model in Figure 10b is a scaled (by 0.95) and smoothed version of the horizontal average of the true model. In this example, except for the maximum offset (3 km), other survey parameters are the same as in the previous example. Figure 11 displays two common-shot gathers simulated in the true model with a 15-Hz Ricker wavelet as the source, and only reflection data will be used for inversion.

Figure 12a is the RTM image migrated with the true model; accordingly, Figure 12b is the migration image for the initial model. Because the initial model significantly deviates from the true model, all the reflectors in Figure 12b, except for the water bottom, are located at inaccurate positions, and some deep reflectors are misplaced by more than the dominant wavelength. In addition, two major faults and two diffractors are poorly imaged in the RTM with the incorrect smooth initial model.

In this example, conventional FWI can only add some details in the scale of one wavelength or less, as indicated by the updated model shown in Figure 13a, which is the inversion result after 10 iterations of conventional FWI. Because of the lack of low wavenumbers in the model (Figure 13a), the corresponding RTM image shown in Figure 13b does not provide any noticeable improvement over the initial result in Figure 12b. Figure 14a and 14b show two shot gathers simulated in the updated model shown in Figure 13a. Compared with the true gathers in Figure 11, we observe that the synthetic data match the true data well, but only for small offsets. As illustrated in Figure 14c and 14d, traveltime shifts estimated by dynamic warping between the true data and the synthetic data indicate large misalignments at far offsets. Without low-wavenumber updates, conventional FWI cannot avoid cycle skipping, even after 10 iterations.

However, after one iteration FWI adds reflectors to the initial model at possibly incorrect positions, WERTI begins to update the low-wavenumber background model and to correct the reflector positions. These two parts take place in the same iteration (step 1) but in an alternating fashion. To ensure that WERTI updates only the background model, we apply a spatial low-pass filter to the gradient

Figure 15a shows the updated low-wavenumber background model after 10 iterations of WERTI. With this background model, we can significantly improve the RTM image quality, as illustrated by Figure 15b. Compared to the initial migration image in Figure 12b, Figure 15b shows that two major faults and diffractors are better resolved, and reflectors are now close to the true positions shown in Figure 12a. We plot the convergence of WERTI in this example in Figure 16. The traveltime misfit function (Figure 16a) decreases rapidly to nearly zero after WERTI updates the low-wavenumber background and corrects reflector positions. Because this is a synthetic example, we are able to also plot a model misfit curve (Figure 16b), which also decreases significantly with low-wavenumber updates.

In step 2, we can further improve the WERTI result in Figure 15a by using it as a new initial model and begin iterations of conventional FWI. In this way, we can refine the low-wavenumber WERTI model with high-wavenumber details. Figure 17 shows the final model after 10 iterations of FWI. With the addition of high-wavenumber details, the final model represents most of the major features including faults and diffractors. We once again simulate synthetic data in the final model and show two gathers in Figure 18a and 18b. Compared with the true data in Figure 11a and 11b, seismograms are well aligned at the near and far offsets. The small traveltime shifts displayed in Figure 18c and 18d illustrate this alignment in a quantitatively compelling way. Finally, we compare the convergence of FWI after WERTI to conventional FWI in Figure 19. Conventional FWI can decrease the data misfit function, but due to cycle skipping, the model misfit function in fact does not decrease. On the other hand, FWI after WERTI decreases the data misfit and model misfit curves more rapidly.

## DISCUSSION

Conventional FWI with reflection data requires a good initial model to avoid cycle skipping. In contrast, the initial model in WERTI can be far away from the true one, so that cycle skipping occurs in high-frequency reflection data. As shown in previous examples, WERTI with dynamic warping can minimize traveltime shifts between synthetic and recorded reflections, even though they initially deviate from each other by more than half a cycle.

The dynamic warping optimization problem implied by equation 7 is similar to the FWI problem in equation 12. Dynamic warping finds globally optimal time shifts, instead of velocities in FWI, to solve a nonlinear least-squares problem in the objective function (equation 7). Because dynamic warping is a global optimization, it is less sensitive to noise and phase errors in wavelet estimation than correlation-based methods. Figure 2 and more examples shown by Hale (2013) demonstrate the robustness of DIW when used to estimate time shifts from data with other differences unrelated to time shifts, such as noise and waveform differences. Still, further testing is necessary to assess the ability of DIW to estimate time shifts between synthetic and recorded data in the context of WERTI, especially when WERTI is applied to real data. In this paper, WERTI minimizes only time shifts between synthetic and recorded data. Future extension of WERTI to minimize spatial shifts along the source and receiver coordinates of data may also be necessary.

Each WERTI iteration is more expensive than each FWI iteration; the additional cost is primarily due to the new gradient

## CONCLUSION

In this paper, we proposed WERTI to estimate the low-wavenumber components of subsurface velocity models using seismic reflections. WERTI requires no picking, and it uses dynamic warping to detect time-varying time shifts, without suffering from cycle skipping problems that often plague correlation-based methods. By minimizing only traveltime errors, WERTI can overcome local minima that cannot be avoided in reflection FWI with insufficient low frequencies in recorded data.

We design an inversion process by combining WERTI and FWI to estimate the low- and high-wavenumber components of the velocity model alternately. Tests on two synthetic models demonstrate that this process can effectively recover subsurface velocity models using only high-frequency reflection data. This alternation strategy can be extended to invert for the low-wavenumber background velocity and the high-wavenumber impedance contrast, which in fact characterize the subsurface more accurately and, of course, requires density in the state and adjoint-state wave equations.

## ACKNOWLEDGMENTS

We are grateful to Y. Luo (Saudi Aramco) for insightful discussions and suggestions. We thank the associate editor S. Operto, R.-É. Plessix, G. Lambaré, and one anonymous reviewer for their thorough review and many constructive suggestions that led to significant improvements of the paper. This work is sponsored by a research agreement between ConocoPhillips and the Colorado School of Mines (SST-20090254-SRA).

## APPENDIX A ADJOINT-STATE METHOD FOR WERTI

In this appendix, we follow the adjoint-state approach described by Plessix (2006) to derive the gradients

Recall that DIW aims to find

Equations 2, 3, and A-1 are the three state equations in WERTI. We define an augmented functional

Adjoint-state variables

Finally, we can obtain the gradients in equations 10 and 11:

- Received January 4, 2013.
- Revision received June 13, 2013.
- Accepted June 30, 2013.