We present an algorithm which allows us to model wavefields with frequency-domain methods using a much smaller number of frequencies than that typically required by the classical sampling theory in order to obtain an alias-free result. The foundation of the algorithm is the recent results on the compressed sensing, which state that data can be successfully recovered from an incomplete measurement if the data is sufficiently sparse. Results from numerical experiment show that only 30% of the total frequency spectrum is need to capture the full wavefield information when working in the hard 2D synthetic Marmousi model.
Seismic wavefield modeling is a technique which can be carried out in either the time domain or the frequency domain. When solutions to the Helmholtz equation for all frequencies are available, modeling in the frequency domain is equal to modeling with time-stepping methods (Symes (2007)). However, if we fail to accurately obtain solutions for the entire discretized spectrum of the signal, aliasing artifacts will be visible in the resulting modeled wavefield.
We are motivated by the long-standing observation that seismic imaging algorithms formulated in the frequency domain can
yield useful alias-free images by acting only on a very limited number of frequencies; see Sergue and Pratt (2004) and Mulder and Plessix (2004). In this respect, seismic inversion in the frequency domain is considered more serviceable compared to time-stepping, where by definition every time step must be considered in order to obtain results. Sergue and Pratt (2004) showed that the extent to which we can exploit this property of frequency-based inversion is related to the range of available offsets in the seismic data, and furthermore presented a strategy for selecting an optimal subset of frequencies to consider.
Following the spirit of frequency-based seismic inversion, we believe that a similar property of information redundancy exists for the seismic wavefield itself in the frequency domain, which we can apply to save computation in seismic wavefield modeling. As we shall see, our conjecture follows naturally from considering a subset of monochromatic Helmholtz equation solutions as a restricted sampling of the modeled seismic wavefield in the Fourier domain, which has been demonstrated earlier by Lin and Herrmann (2007). In this approach, the mathematical machinery of compressed sensing ( Cand`es et al. (2006); Donoho (2006)) naturally leads to a robust method of anti-aliasing the modeled wavefield in the time domain by solving a de-noising problem, an idea introduced in Hennenfent and Herrmann (2006).
In this paper, we describe an algorithm for computing acoustic wavefields by using only a limited number of frequencies. The algorithm is based on sparsity-promoting recovery formulated in the context of compressed sensing. We also present experimental results showing that this algorithm can compute wavefield using only 30% of the total frequency information even in models with sharp discontinuities, such as the hard 2D Marmousi model.
Solving the Helmholtz equation is the standard way to model the propagation of seismic wavefield in the frequency domain through an inhomogeneous media.