An Efficient Approach for Computing the Frequency Response of Seismic Waves in Heterogeneous, Anisotropic Viscoelastic Media With FDTD+PSD Modeling

Computation of the frequency response (phase and magnitude) of seismic waves propagating in heterogeneous, anisotropic, viscoelastic media is required for a number of seismology efforts, including frequency-domain full-waveform inversion and basin modeling. When the frequency response is required at a limited number of locations, it can be computed efficiently with a memory variable, staggered grid finite difference time domain (MV-SG-FDTD) code by taking the Fourier transform of the received waveforms generated by a broadband source. However, when the frequency response is required at many or all grid locations, as in frequency-domain full-waveform inversion, the memory requirements for storing the waveforms make this approach prohibitive. An alternative approach is to compute the frequency response by formulating the SG-FDTD equations in the frequency domain. For O(4) accuracy spatial differencing, the resulting 2-D system of implicit equations for the unknown particle velocities and stresses at the finite difference cell locations yields a complex system matrix of order equal to [5xMxN]², where M and N are the number of finite difference cells in the x and z directions, respectively, and the 5 comes from the five unknowns in a 2-D staggered grid finite difference stencil (2 particle velocities, and 3 stresses). The complex matrix equation has the following properties: (1) sparse (total number of unknowns equals [45xMxN-20xM-20xN]), (2) banded (4 bands of sparse 5x5 submatrices on each side of the main diagonal), and (3) non-Hermetian. We found that the resulting system matrix is typically much too large for direct matrix solvers such as sparse LUD for seismic problems of reasonable size. Additionally, our attempts to solve the frequency domain system of implicit equations using iterative sparse matrix solvers, such as QMR (quasi minimum residual method), tended to suffer from poor convergence problems when using standard preconditioning (e.g., Jacobi). In this presentation, we demonstrate that an efficient approach for computing the frequency response of a heterogeneous, anisotropic, viscoelastic medium can be achieved by running an explicit MV-SG-FDTD code with a harmonic wave source out to steady-state, and then extracting the magnitude and phase from the transient data via a phase sensitive detection algorithm (PSD). The PSD algorithm requires integration over only several cycles of the waveform to obtain accurate phase and magnitude estimates. Because this integration is performed on-the-fly, there is no need to store waveforms at the grid locations. Preliminary tests indicate that it should be possible to superimpose multiple frequencies at the source and extract the magnitude and phase for each frequency using the PSD algorithm, opening up the potential for obtaining the multi-frequency response of a heterogeneous, anisotropic, viscoelastic medium with a single FDTD run.