My main professional interest is seismic imaging of the earth interior by tomography and full waveform inversion (FWI).
I am interested in theoretical, algorithmic and application aspects of FWI. FWI is a PDE-constrained optimization method which aims to estimate mechanical properties of the earth interior by fitting the waveforms of seismic records (Figure 1).
Figure 1: Illustration of the data fitting procedure. The recorded seismograms are in black, the simulated seismograms are in red. From left to right, the simulated seismograms are computed in a subsurface model obtained by traveltime tomography, mono-parameter FWI and multi-parameter FWI accounting for density and attenuation. Waveform matching improves from left to right as more complete information in the seismic data is exploited.
The PDE constraint is the wave equation where the mechanical properties of the earth are embedded in its spatially-varying coefficients. In its conventional form the objective to be minimized is the least-squares norm of the difference between the seismic observables and the restriction of the numerically-simulated wavefield at the receiver positions. We classically solve this optimization problem with local optimization techniques (gradient-based methods) and a reduced-space (variable projection) approach: we transform the PDE-constrained optimization problem into an unconstrained problem by enforcing the exact solution of the wave equation for the current guess of the subsurface parameters in the objective and solve the unconstrained problem with Newton-type algorithm where the gradient of the objective is computed with the so-called matrix-free adjoint-state method. It is well acknowledged that this inverse problem is highly nonlinear due to the oscillatory nature of seismic waves. The linearization of the FWI relies on the single-scattering Born approximation which requires the initial subsurface model to allow for the fitting of the recorded traveltimes with an error that does not exceed half the period. This kinematic-accuracy condition is challenging to achieve for real-life problems in particular for long-offset acquisitions (Figure 2).
Figure 2: Sensitivity kernel of FWI in the frequency domain for one source-receiver pair, one frequency and a homogeneous background model. The interference between the source and receiver wavefield generates the "iso-phase" surfaces, also called isocrhones associated with single-scattered paths connecting the source to the receiver via scattering points in the subsurface. These iso-phase surfaces are in fact Fresnel zones, defined by a subset of positions x in the subsurface such that the paths S-x-R have traveltimes that do not differ by more than half a period. In FWI, the residuals (the difference between the recorded and the simulated seismograms) are backprojected onto these Fresnel zones provided that the modeled and the recorded components of the residual have traveltimes that do not differ more than half a period (if this condition is not fulfilled they will be indeed backprojected unconsistently onto two different isochrones and FWI will converge to a poor minimizer. In the figure, k denotes the wavenumber component locall mapped in the subsurface at the position of the scatterer. The modulus of this wavenumber vector corresponds to the width of the isocrhone and is related to the local wavelength and the cosine of half the scattering angle. This relationship has been established in the framework of diffraction tomography. Note how the width of these isocrhones decrease as the scattering angle decreases.
My current interestes are to overcome this pitfalls with two different approaches:
The first aims to improve the accuracy of the initial subsurface model that is built by traveltime + slope tomography. Slope tomography seeks to track in a semi-automatic way locally-coherent events in the seismic data volume to perform a dense picking of traveltimes and slopes (the horizontal component of the slowness vectors ideally at the source and receiver positions). This dense picking is amenable to high-resolution velocity model building. Traveltimes are sensitivie to the long-wavelength distribution of the velocity field, while the slopes carry out a complementary information on the velocity gradients. We are currently working on a new formulation of slope tomography based on eikonal solver and the adjoint-state method. We aim to apply this technology to towed-streamer data as well as to multi-component ocean-bottom seismometer (OBS) data.
These developments have been performed during the PhD of Borhan Tavakoli and continue during the ongoing PhD of Serge Sambolian.
The second approach seeks to extend the linear regime of the FWI by extending its search space. Said in another way, we develop approaches to foster the data fitting and prevent the cycle skipping pathology accordingly by relaxing the wave equation constraint (Figure 3). We implement this idea by implementing FWI with the alternating direction method of multipliers (ADMM) i.e., an augmented Lagrangian method with operator splitting. This method is implemented with Total Variation regularization. This method was originally developed by T. van Leeuwen and F. Herrmann with penalty method and H. Aghamiry has recently improved the wavefield reconstruction method by replacing a penalty method by the alternating direction method of multiplier (which is equivalent to the Bregman method). This makes the optimization method far less sensitive to the penalty parameter and improves dramatically the convergence.
These developments are performed during the PhD of Hossein Aghamiry under the co-supervizion of Professor Ali Gholami (University of Tehran) and myself.
Figure 3: Illustration of the wavefield-reconstruction method. The true subsurface model is shown on the top left panel. It shows a rectangular inclusion in a velocity-gradient background. Let's imaging that we start our imaging procedure from the homogeneous velocity model shown on the bottom left panel. The top right panel shows the difference between the wavefield computed in the true subsurface model and the initial counterpart, that is the wavefield scattered by the difference between the two subsurface models shown in the left. The middle right panel shows this scattered wavefield predicted by the initial model (indeed, it is zero). The bottom right panel shows the scattered wavefield when we reconstruct the wavefield that best satisfies in the least-squares sense the wave equation solved in the initial model and fits the data. The scattered wavefield is now closer to the scattered wavefield shown on the top panel. If we update the subsurface parameters from this wavefield residual, we extend the linear regime of the FWI (Figure from Hossein Aghamiry).
From the algorithmic viewpoint, we are developing at Geoazur in-house slope tomography code as well as 3D time-domain and frequency-domain FWI codes. The time-domain FWI code is developed by L. Combe and the frequency domain code by A. Miniussi and myself. In this latter code, we use the MUMPS multifrontal solver to solve efficiently the time-harmonic wave equation with multiple right-hand sides.
We are working on different applications on OBC data from the North sea (Figure 4 and Publications 4, 12, 13 and 17) and on OBS data from the eastern Nankai trough in Japan (Publication 8) in collaboration with Andrzej Gorszczyk and Michal Malinowski (Institute of Geophysics of the Polish Academy of Science).
Figure 4: Velocity model of a North Sea oil field obtained by 3D frequency-domain FWI (Operto et al., 2015).
I am also interested in application of FWI on earthquake data, in particular on teleseismic data.
Stephen Beller and Vadim Monteiller have developed a 3D FWI code for lithospheric imaging from teleseismic data. A sensitivity analysis of the data to the acquisition geometry, the subsurface parametrization and the initial model was performed as well as an application on the data from the CIFALPS experiment in the western Alps (publications 6 and 7).
I am launchind the 3DWIND project which aims to develop the FWI technology for deep crustal imaging from multi-component OBS data. This raises three important challenges related to the acquisition design of OBS survey amenable to slope tomography and FWI, high-performance computing challenges to address FWI applications involving up to billions of parameters and methodological challenges related to the nonlinearity of the FWI for long-offset deep-crustal stationary recording geometries.
I am also involved in the SEFASILS project which aims to perform a regional seismic exploration of the Northern Ligurian margin in November-December 2018.
First-arrival slope tomography (PhD thesis of Serge Sambolian)
