An Iterative Multilevel Method For Computing Wavefields In Frequency-domain Seismic Inversion

Erlangga, Yogi A. (University of British Columbia) | Herrmann, Felix J. (University of British Columbia)



We describe an iterative multilevel method for solving linear systems representing forward modeling and back propagation of wavefields in frequency-domain seismic inversions. The workhorse of the method is the so-called multilevel Krylov method, applied to a multigrid-preconditioned linear system, and is called multilevel Krylov-multigrid (MKMG) method. Numerical experiments are presented for 2D Marmousi synthetic model for a range of frequencies. The convergence of the method is fast, and depends only mildly on frequency. The method can be considered as the first viable alternative to LU factorization, which is practically prohibitive for 3D seismic inversions.


Imaging of the earth’s subsurface can computationally be done in the frequency domain by inverting seismic data using the adjoint-state method; see, e.g., Plessix (2006). One advantage, among others, of doing inversion in the frequency domain is the possibility of inversion with only a small subset of random frequencies, as shown, e.g., in Sirgue and Pratt (2004) and Mulder and Plessix (2004). In relation to this, it is demonstrated by Lin et al. (2008) (in this proceedings) that it is possible to recover the wavefield with only 30 per cent of frequencies required by the classical Nyquist sampling theory, by solving a sparsity-promoting recovery problem in, e.g., the Fourier domain. In combination with less number of grid points required to resolve the wavefield, seismic inversion in the frequency domain will now require less complexity as compared to the time domain migration.

Translated to the seismic language, the adjoint-state method consists of forward modeling, in which a wave equation is solved for a given velocity background and source position, and computing a correction to the given velocity background. The correction is determined by minimizing the misfit between the seismic data, recorded at the receiver positions, and the forward model data at the positions corresponding to the physical receivers. The gradient of the misfit functional is closely related to the migrated image in geophysical applications, and can be computed by multiplication of forward and back propagated wavefields. These wavefields are computed via a finitedifference approximation of the wave equation. As mentioned, e.g., in Pratt et al. (1998), wavefield computations may be done by direct methods, by first constructing the LU factors of the finite-difference wave equation matrix A.

Iterative methods are usually considered as an alternative to direct methods (Saad (2003)). An iterative method relies mostly on matrix-vector multiplications with A. Hence, an iterative method is less demanding in terms of memory and computational complexity. Iterative methods are, however, generally less robust as compared to direct methods, and in particular for wave simulations, all efficient, ad hoc iterative methods fail to converge for frequency as high as 5 Hz, in the case of, e.g., the 2D Marmousi synthetic model. Erlangga et al. (2006) proposed an iterative method based on Krylov subspace method, preconditioned by a damped acoustic wave operator. The preconditioner is inverted approximately by one multigrid iteration. In this paper, we call this method the MG method. While the method improves the convergence significantly, the convergence still shows dependence on angular frequency.