Finite-difference modelling is a common method of solving wave equations and is the foundation for several of Acceleware’s products. One problem with finite-difference methods is that material interfaces that do not align with the rectangular grid are not easily represented. At the very least, naïve implementations may represent a structure different from the one intended. The typical example is a sloping interface between two different materials. If the grid is sampled based on nearest neighbor interpolation, then the model will have a staircased slope, as shown in Figure 1.
Figure 1 - Staircase representation of a flat sloping interface
In our electromagnetic products, a combination of sub-gridding and material averaging is used to accurately represent man-made structures, but in the seismic industry, extremely sparse grids are used to reduce the memory required for modelling. Therefore, staircasing errors are a common problem.
One solution I was able to apply originates from Rune Mittet’s work titled “On the Internal Interfaces in Finite Difference Schemes”  and some internal conversations with our Chief Science Office, Michal Okoniewski. The key is to keep in mind the continuous representation of a discrete set of samples. In the case of a sharp material interface, simply jumping from one value to the next results in oscillations known as the Gibbs phenomenon. Instead, we need to represent the jump accurately within the available bandwidth (in the wavenumber domain).
In realistic cases, there is more than one sharp interface, interfaces intersect one another in complicated ways, etc. Therefore, we assume only that the material parameters are known to much greater accuracy than the computational grid. In that case, a rectangular grid much finer than the computational grid can be generated. It is then filtered in the wavenumber domain such that it can be safely decimated to the computational grid. Doing this, results in substantially less staircasing errors, as shown in Figure 2 and 3.
Figure 2 - Synthetic modelling of a 2D model comprised of sloping reflectors. When the model is filtered to the appropriate bandwidth, the staircasing errors are significantly reduced. Multiple shots similar to this result were used to generate the images in Figure 3.
Figure 3 - RTM of the modelled result superimposed on the density model. The upper image uses a staircased model whereas the bottom one is filtered to the appropriate bandwidth. Note the middle reflector especially is represented as a series of flat reflectors at different depths instead of a smooth slope. The circled sections are what can be viewed below in Figure 4.
Figure 4 - Zoomed in sections from Figure 3, demonstrating the difference between a staircased model and one which is filtered to the appropriate bandwidth.
Unfortunately, filtering using discrete Fourier transforms uses an impractical amount of memory because computational grids already fill a large percentage of available memory. Additionally, sampling 10 times finer in 3D requires 1000 times as much memory. Instead of filtering in the wavenumber domain, I applied a filter on a finely sampled grid in the spatial domain. Since the desired output is the coarsely sampled computational grid, the finely sampled grid does not need to be kept in memory. Instead, a sliding window of samples are evaluated, filtered, and then decimated. This approach is easily multithreaded, and only requires redundant evaluations where the finely sampled grid is partitioned over threads.
In summary, staircasing errors can be significantly reduced by sampling model parameters that have been filtered to the appropriate bandwidth. In order to represent arbitrary models, a finer, but staircased, model is generated, filtered, and then decimated. These 3 steps are combined in order to make the process computationally feasible.