An Helmholtz iterative solver for the three-dimensional seismic imaging problem?

More Info


During the last three decades high-frequency approximations or paraxial (one-way) approximations of the wave equation have been successfully used to process seismic data in two-dimensional or three-dimensional spaces. With the increase of computer power and the need to take into account complex geological structures, solutions of the exact wave equation either in time domain or in frequency domain are required. The time-domain wave equation can be solved with a marching approach in time and an explicit scheme. Currently three-dimensional imaging algorithms, so-called reverse time migrations, are developed based on time-domain solutions. The frequency-domain wave equation, namely the Helmholtz equation in the acoustic case with constant density, leads to a large sparse linear solver. In a two-dimensional space a direct solver, based on a LU decomposition, can be used. This leads to attractive and efficient imaging methods faster than the time-domain imaging methods. In a three-dimensional space, due to the size of the matrix, an iterative solver is required. In practice three-dimensional frequency-domain solutions are not used because most of the iterative solvers are not robust at seismic frequency or are too expensive. The difficulty arises because the matrix of the linear system is not definite and because the model size generally corresponds to several hundreds of wavelength in the three space directions. A new preconditioner based on a heavily damped wave equation and a multi-grid solver has been proposed by Erlangga et al.. The iterative solver is a preconditioned bicgstab. The preconditioner corresponds to the approximated inverse of the heavily damped wave-equation after one multigrid cycle. Numerical results have shown that this new iterative solver is robust. However does a frequency-domain imaging algorithm based on this solver compete with a time-domain imaging algorithm? To handle realistic three-dimensional geophysical problem the multigrid solver is implemented with a standard full-weighting coarsening, a tri-linear interpolation, a matrix-free approach, and an 8th order scheme in space. Three-dimensional numerical results indicate that the number of bicgstab iterations linearly increase with frequency. Based on this result, the number of the floating point operations needed for the time-domain and the frequency-domain imaging algorithms is roughly estimated. It is then shown that the migration algorithm is about one order of magnitude faster in the time domain and that the velocity modeling building algorithm based on a least-squares formulation and a scale separation is about three time faster in the frequency domain.