Exact local refinement using Fourier interpolation for nonuniform-grid modeling

1.   Introduction

Numerical simulation is necessary to understand complicated wave phenomena and to further estimate the physical properties (e.g., inner geometrical structures or even velocity distributions) of the media with the aid of a computer. The finite-element, spectral-element and grid method (Zhang JF and Gao HW, 2009) can freely handle arbitrary geometry structures using a globally varying grid. Whereas, their resolution is low when using a coarse grid; on the other hand, the computational cost would be too large with a fine grid since it involves calculating the inverse of a large-scale matrix. In contrast, the numerical solver using a uniform grid (e.g., the finite-difference method and pseudo-spectral method) has less memory demand and computational cost, since it does not need to store the coordinates of the grid system and the wavefield can be explicitly updated without solving the inverse of the matrix. Thus, the numerical solver in the time domain using a uniform grid is one of the most important methods, and it is popular in various practical applications due to its simplicity and good performance (Etgen and O’Brien, 2007; Liu Y and Sen, 2013; Zhang JH and Yao ZX, 2013; Tan SR and Huang LJ, 2014).

However, the numerical solver using a uniform grid will encounter serious problems in the presence of tiny structures in large-scale media. For example, the influence of tiny structures on the wave propagation is fairly sensitive to the thickness of the tiny structure (Zhang JF and Gao HW, 2009), as shown in Figure 1. If we divide the model by a uniform grid interval (e.g., 0.1 m) that is consistent with the minimum thickness of the tiny structures, we can only handle small-scale models; for a large-scale model, the memory demand and computational cost would be too large for most current computers, especially in 3D cases. For example, the total grid number of a 1 km³ cubic is up to 10¹² when using a uniform grid interval of 0.1 m.

Figure  1.  Velocity model with tiny structures (a) and its impulse response (b). The velocity of the homogenous background is 3000 m/s, and the velocity within the tiny structures is 1000 m/s. The vertical tiny structure ranges from (45, 40) to (45.1, 60); and the horizontal tiny structure ranges from (46.2, 85) to (51.2, 85.2). The black star shows the position of the source. The source is a Ricker wavelet with a dominant frequency of 1000 Hz. The grid interval is 0.1 m and the temporal interval is 0.01 ms. We solve the scalar wave equation with the pseudo-spectral method. The snapshot of the wavefield shown in (b) is at 12 ms.

DownLoad: Full-Size Img PowerPoint

A reasonable way to reduce both memory demand and computational cost is to use a fine grid over target zones but still use a coarse grid over other zones (Moczo, 1989; Jastram and Behle, 1992; Falk et al., 1998; Tessmer, 2000; Sergio and Oliveira, 2003; Zhao HB and Wang XM, 2008; Song GJ et al., 2010; Chu CL and Stoffa, 2012; Zhang ZG et al., 2013). However, due to using a nonuniform grid, these methods rely on either wavefield interpolation (i.e. upsampling) or varying weighting coefficients using Taylor expansion (Huang C and Dong LG, 2009). We know that Taylor expansion has high accuracy only when the variables are close to the expansion point; thus, the finite-difference method using Taylor expansion would encounter numerical dispersions when using a nonuniform grid within a thick transition zone between a coarse grid and a fine grid. In contrast, the wavefield interpolation is more popularly used to upsample wavefields from coarse grids to fine grids (Falk et al., 1998; Tessmer, 2000; Song GL et al., 2010; Zhang ZG et al., 2013). Unfortunately, most current wavefield interpolation methods (e.g., Lagrange interpolation, spline interpolation, or Hermit interpolation) are not accurate enough and, thus, would introduce apparent artifacts around the transition zone. Generally, the wavefield interpolation can only handle a small upsampling ratio (e.g., usually 2–5), and would fail when the upsampling ratio is bigger than 10 due to too abrupt wavefield perturbations. We need to develop a completely new kind of wavefield interpolation that is accurate enough for a high upsampling ratio. Recently,Kostin et al. (2015) proposed a novel approach to local time-space grid refinement, which is able to handle a high upsampling ratio since the local refinement of the wavefield on a coarse grid is based on the Fourier interpolation. Whereas, their method is valid only for an odd upsampling ratio. In addition, both the theoretical basis and implementation details of Fourier interpolation are still not concise and clear in the literature.

In this paper, we propose a new scheme for upsampling the wavefield using Fourier interpolation. Figure 2 illustrates the flowchart of Fourier interpolation, and Figure 3 shows the procedure for updating the wavefield between a coarse grid and a fine grid. A significant advantage of the Fourier interpolation is that it would be theoretically exact if the input wavefield is band limited, since it reconstructs temporal samples by padding zeros to the spectra of the wavefield. On the other hand, the Fourier interpolation does not have a limitation on the upsampling ratio, which is fairly attractive for our purpose. In addition, this technique is easy to implement since it only involves three steps: first, a forward Fourier transform; then, inserting zeros at the Nyquist frequency; and finally, an inverse Fourier transform.

Of course, there are some aspects that should be noted to guarantee both stability and accuracy. First, we need to apply a smooth enough window function or a taper zone outside the wavefield that we hope to refine, which is to make sure the fetched wavefield is periodic after tapering the boundary. A smooth enough window function usually means involving more points around the target zone. Second, the shape of the smooth window function should be rectangular (cuboid in 3D) to facilitate a fast Fourier transform. Third, when the source or receiver is close to the target zone, this scheme would encounter problems since the amplitude of the wavefield would be tapered by the smooth window function; fortunately, we can completely enclose the source or receiver by enlarging the window function, although this would apparently increase the computational burden.

Figure  2.  Flowchart of the Fourier interpolation. (a) 1D case (i.e. the refinement of sampling); (b) 2D case (i.e. the enlargement of a picture). (a1) The original data that are ready to be refined or interpolated, whose total number of samples is N; (a2) the amplitude spectrum of the original data; (a3) the procedure of splitting the Fourier transformed data at the Nyquist frequency f_N and inserting (M–1) times zeros, where M×N is the total number of samples after interpolation; (a4) the interpolated data (i.e., the real parts of the inverse Fourier transform after inserting zeros). Circles and dots in (a4) are original and interpolated samples, respectively. We use M=3 for (a) and M=2 for (b), respectively.

DownLoad: Full-Size Img PowerPoint

Figure 3 illustrates the procedure of updating the wavefield between a coarse grid and a fine grid. First, we only consider the wavefields in the background by intentionally omitting the tiny fractures in the target zone, which can avoid artifacts due to too poor spatial sampling around tiny structures on the coarse grid. Second, we perform regular forward modeling on the coarse grid but for only one time step. Third, we fetch the wavefield on the target zone with a smooth window function at both current and previous time steps. The wavefield fetched from the coarse grid would be first interpolated and then assigned to the corresponding updating points on the fine grid. Note that the wavefield exchange happens only on the updating points shown in Figure 3. Then, we perform regular forward modeling on the fine grid for the target zone until the time step is consistent with that used by the coarse grid. Next, we fetch the wavefield on the fine grid with a smooth window at the current time step. Finally, the wavefield fetched from the fine grid would be assigned to the corresponding updating points on the coarse grid, as shown in Figure 3. We should use a wide enough taper zone (e.g. 50–100 points) for the smooth window to reduce possible artificial reflections; meanwhile, the taper zone should be far enough (say 10 points) away from the target zone.

Figure  3.  Illustration of local refinement using Fourier interpolation. (a) for 1D case; (b) for 2D case. The circles are on the coarse grid, which are around the fine grid. The wavefields on the coarse grid are extracted from the total wavefield at the current and previous steps, respectively. The local fine grid should cover the target zone and a proper taper zone is necessary to avoid artificial boundary reflection. The coarse grid should cover the target zone and also requires a taper zone to avoid artifacts.

DownLoad: Full-Size Img PowerPoint

2.   Methodology

For an input time series $x\left[ n \right),\;n = 0,1, \cdots ,N - 1,$ its Fourier transform is defined as (Sacchi et al., 1998)

The inverse Fourier transform is defined as

In equations (1) and (2), ${\rm i} = \sqrt { - 1}$ , N is the total number of samples, n is the sample index in the temporal domain, and k is the sample index in the frequency domain. If the uniform interval between the samples of $x\left[ n \right]$ is $\Delta t,$ the total length of the time period is $T = N\Delta t.$ According to the sampling theory (Oppenheim et al., 1999), the maximum non-zero frequency component of a band limited signal would be smaller than the Nyquist frequency to avoid aliasing, where the Nyquist frequency is defined as

In other words, the amplitude spectrum of a band limited signal must be zero around the Nyquist frequency. Thus, if we pad $\left( {M - 1} \right) \times N$ zeros at the Nyquist frequency, we could not add any new frequency component into the spectrum and could not remove any existing components from the spectrum; the only change is that the total number of samples becomes $M \times N$ , or the Nyquist frequency is extended to a new position ${f'_N} = M{f_N}$ . Therefore, we have a new time interval

This indicates that padding zeros in the frequency domain at the Nyquist frequency is equal to reducing the temporal interval, which is basically what an interpolation does. Different from traditional interpolation methods (e.g., polynomial interpolation, and spline interpolation), this procedure is exact since it does not rely on any assumption except that the signal should be band limited and sampled with a uniform interval. Note that the inverse Fourier transform becomes

The flowchart of the Fourier interpolation is illustrated using both 1D and 2D data, as shown in Figure 2.

3.   Experiments

The 2D scalar wave equation reads

where t is the time axis, $v\left( {x,z} \right)$ is the velocity of wave propagation, and $u \equiv u\left( {x,z;t} \right)$ is the wavefield. We use the Ricker wavelet

as the source time function, where $\omega$ is the dominant frequency. We solve the wave equation using the pseudo-spectral method (Kosloff and Baysal, 1982) to illustrate the local refinement of the wavefield using the Fourier interpolation. The wavefield at the next step is iteratively predicted by the wavefield at the current and previous steps as follows:

where $F_{\! x}^ +$ and $F_{\! z}^ +$ are 1D forward Fourier transforms along the x- and z-directions, respectively; $F_{\! x}^ -$ and $F_{\! z}^ -$ are 1D inverse Fourier transforms along the x- and z-directions, respectively; ${k_x}$ and ${k_z}$ are the wavenumbers along the x- and z-directions, respectively.

The wavefield on the coarse grid is extracted from the previous and current time steps, respectively, using a tapered window function (e.g., Hanning window with 100 points) over the target zone, as shown in Figure 3. Then, the extracted wavefield is refined by the Fourier interpolation shown in Figure 2. The refined wavefield is assigned to the local fine grid, and we can independently perform a local forward modeling on the local fine grid using Δt/M. Such a temporal interval for a local fine grid can guarantee the stability conditions and can coincide with the temporal interval used by the coarse grid. When the current time coincide with the time for the coarse grid, we update the wavefield on the coarse grid by extracting the wavefield at the same positions that are generated on the local fine grid.

To verify the accuracy of the proposed local refinement using Fourier interpolation, we compare our method with some typical interpolation methods. We perform a forward modeling using the pseudo-spectral method with a grid interval of 0.01 m. The wavefield at 0.001 ms is taken as a reference of fine grid data, whose sample number is 2000×2000. We extract one sample every 50 points along the x- and z-directions (i.e., M=50), respectively. We use these fairly sparse data as the input of the Fourier interpolation, whose sample number is 40×40. Obviously, traditional interpolation methods listed are not able to accurately reconstruct the waveforms especially on the waveform peaks, as shown in Figure 4. The resulting visible waveform errors would lead to apparent discrepancy between the wavefield generated on a coarse grid and that on a fine grid, which would be represented as strong artifacts around the transition zone. In contrast, the biggest error for our method is about two orders of magnitude smaller than those of traditional interpolation methods. Note that, the total number of input samples is only 40×40, which is almost close to the lower limit for the Fourier interpolation. If we have a much bigger number of input points (say 100×100) using either a much wider window function or a much smaller spatial interval, the accuracy can be further improved.

Figure  4.  Comparison of local refinements between different methods. (a) zoom in of refined wavefield; (b) waveforms; (c) waveform errors. There are five parts in (a) that are separated by solid lines. The fine data (the lower part of the left-upper quadrant) are used to check the accuracy of different interpolation methods. The bold dashed segment in the right-bottom quadrant of (a) is the profile where the waveforms shown in (b) are extracted from.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

The strictest requirement for the Fourier interpolation is that the fetched wavefield is band limited. We know that the number of sampling points in each wavelength is usually at least several for most finite-difference methods to avoid numerical dispersions. Even for the pseudo-spectral method and the finite-difference method proposed by Yang DH et al. (2003), which are almost free of numerical dispersion, the number of sampling points in each wavelength for the highest frequency component is still bigger than 2, since 2 is the lower limit for all methods. This indicates that the wavefield generated by numerical simulations is smooth enough spatially, which has no abrupt jump or discontinuity; in other words, the wavefield generated is band limited. Thus, we should only guarantee that the wavefield is periodic after we fetch the data using a smooth window; fortunately, this requirement is easy to fulfill by simply applying enough number of points in the smooth window function. Therefore, the Fourier interpolation is feasible for upsampling the wavefield generated by numerical simulation on a coarse grid, only if we use a proper smooth window. Such a simple requirement would be helpful for the extensive use of this technique.

5.   Conclusions

The upsampling of a wavefield from a coarse grid to a fine grid is essential for numerical simulation methods that use a nonuniform grid. We present a high-accuracy local refinement scheme using the Fourier interpolation. The proposed scheme is superior to traditional interpolation methods since it is theoretically exact if the input wavefield is band limited. Traditional interpolations would fail at a high upsampling ratio while our scheme still works well in the same situation. The high upsampling ratio is important for a successful numerical simulation in the presence of tiny structures and at large scales, especially for the 3D case.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grant No. 41130418) and the National Major Project of China (under grant 2017ZX05008-007). JinHai Zhang thanks the supports from the Youth Innovation Promotion Association CAS (2012054) and Foundation for Excellent Member of the Youth Innovation Promotion Association (2016).