Spatial Filter for the Pseudo-spectral Implementation of Fractional Derivative Wave Equation

The viscoelasticity of the subsurface media varies spatially, and such viscoelasticity can be represented concisely by a wave equation in the form of fractional temporal derivative (FTD). We have developed a strategy for simulating seismic waves propagating through a heterogeneous viscoelastic model. The FTD is transferred to fractional spatial derivatives (FSDs), and the FSDs are implemented through the fast Fourier transform (FFT), for improving the computational efficiency. However, the FFT implementation is not rigorously applicable to the heterogeneous model. In this paper, we have reformulated the FSD wave equation by introducing a spatial-position dependent filter. This spatial filter corrects the errors that are caused by the assumption of non-heterogeneity in the FFT implementation. This formulation appropriately represents the viscoelastic effect in seismic wave propagation, leading to the improvement on the accuracy of numerical simulation.


Introduction
The anelasticity of subsurface media causes seismic dissipation effect including the energy absorption and velocity dispersion, and has a significant impact on field seismic records and subsequent seismic images (Wang, 2008). Conventional models for describing viscoelasticity include Kjartansson's constant-Q model (Kjartansson, 1979), the Kosky model (Kolsky, 1953), the Kelvin-Voigt model, the standard linear solid model (Zener, 1948), the Cole-Cole model (Cole & Cole, 1941) etc. These models can lead to viscoelastic wave equations in the form of fractional temporal derivative (FTD).
The FTD was introduced to describe the viscoelasticity of the Earth media by Caputo (1967). The FTD form was used also for the fractional Zener model to describe the mechanic combination of the viscoelasticity (Liu & Greenhalgh, 2019). However, directly solving the FTD wave equation presents a numerical challenge in seismic wave simulation. The FTD might be solved in the frequency domain, but the computation is extremely intensive because it requires to solve numerous Helmholtz equations (Vasilyeva et al., 2019). An alternative, but still expensive, method is to solve a convolution equation in the time domain (Carcione et al., 2002). The Grünwald-Letnikov expansion can also be used to calculate the FTD, but this expansion requires large computational memory to store the wavefield history, even the expansion is truncated (Podlubny, 1999).
A practical way to reduce the computational cost of the FTD wave equation is that, when the attenuation is weak, one can transfer the FTD to fractional spatial derivatives (FSDs), and then solve FSDs by Fourier pseudo-spectral method, which greatly lowers the computation memory threshold (Carcione, 2010;Chen & Holm, 2004;Treeby & Cox, 2010). The FSDs in the wave equation can also be decoupled into velocity dispersion and energy dissipation, respectively (Zhu & Harris, 2014), and this decoupled wave equation can be further expanded to viscoelastic and tilted transversely isotropic (TTI) media (Qiao et al., 2019;Zhu & Bai, 2019;Zhu & Carcione, 2014). The accuracy of the FSDs for wave simulation in highattenuation media can further be improved by Taylor series expansion (Mu et al., 2021).
When using the pseudo-spectral method to solve the FSDs in the wave equation, it requires that the viscoelastic parameter bðxÞ is a constant in the space.
In practice, the viscoelastic parameter bðxÞ is averaged over the space x, to generate a constant b for the purpose of spatial Fourier transform. This averaging process will cause errors in numerical calculation. Based on the locality principle, Zhu and Harris (2014) propose to interpolate the solutions generated with every single constant parameter for whole model. The spatially varying order of FSD is a function of bðxÞ. Chen et al. (2016) and Xing and Zhu (2019) use either the Taylor expansion or a polynomial approximation to transfer the spatial-varying order FSD into the constant order FSDs, and then implement the pseudo-spectral method directly. Yao et al. (2017) apply the Hermite distributed approximation to transfer the spatial-varying FSD to an integral of the fractional derivative of delta function, and this fractional derivative can be solved locally. But all these schemes aforementioned scarify the efficiency for accuracy to some extent. For example, when using the polynomial approximation to derive the constant-order wave equation (Xing & Zhu, 2019), the resultant wave equation includes six fractional Laplacian operators. This scheme greatly lowers the efficiency, if compared with the original two-operator equation.
In this paper, we propose a strategy to solve the FSDs for wave simulation in heterogeneous media, that would be straightforward in philosophy and simpler in realization. The strategy is to build a spatially varying correction function, and to insert this spatial filter directly into the averaging scheme. Because this spatial filter is frequency-independent, it is efficiently implemented as a coefficient multiplied to the wave equation. Therefore, this filter improves the accuracy and maintains the high efficiency of FFT implementation.

The Equation Formed with Fractional Spatial Derivatives
The wave equation presented in terms of FTD may be transferred to a wave equation formed with FSDs. The ultimate purpose of this paper is to have an efficient implementation of wave equation with FSDs.
We consider the FTD wave equation (Wang, 2016), as the following where u is the scalar wavefield, x 0 is the reference frequency, b is the viscoelastic parameter, and cðbÞ is the viscoelastic velocity, which means the phase velocity at zero frequency. Assuming the medium parameters to be spatially invariant, we can rewrite it in the frequency-wavenumber domain as where k is the wavenumber, andû is the frequencywavenumber domain wavefield. Now, making an approximation based on weak attenuation assumption (Zhu & Harris, 2014) x where c 0 is the reference velocity, Eq.
(2) is approximated to Applying an inverse Fourier transform to Eq. (4), the generalized wave equation in the temporal-spatial domain may be presented as FSDs: where u is the time-space domain wavefield. Equation (6) decouples the dissipation effect, as the C 1 term is from the real part of FTD, which represents the velocity dispersion and the C 2 term is from the imaginary part of the FTD, which represents the amplitude absorption (Zhu & Harris, 2014). In this equation, C 1 ðbÞ and C 2 ðbÞ are assumed to be spatially independent. However, when forming the FSDs for a general viscoelastic case, both terms are spatially variable, because bðxÞ is a spatial function. For the derivation above, the approximation Eq. (3) is a key condition, so that the complex wavenumberk could be replaced with the real wavenumber k ¼k Re . This approximation to the complex wavenumber is made based on the weakattenuation assumption. For the complex wavenumberk, the real and imaginary part may be written analytically as (Wang, 2019) where vðb; xÞ is the phase velocity, aðb; xÞ is the attenuation coefficient, and the absolute value is for satisfying the anti-Hermitian property of the complex wavenumber. The factors A and B were given by Wang (2019). Figure 1 illustrates the accuracy of this weak-attenuation assumption with the ratio R ¼ k=jkj. The reference frequency x 0 is set as 500 Hz, which is the highest possible frequency of exploration seismology (Wang & Guo, 2004), and then x=x 0 1. Figure 1 shows that the ratio R is close to 1, for b 0:75, and the relative error is less than 0.5%. When b 0:351, the ratio R is approaching to 1 with decreasing b. But when b [ 0:351, the ratio R is increasing with an increasing b. Therefore, the weak-attenuation assumption is valid for b 0:351, to the most areas interested by seismic application (Kolsky, 1953;Futterman, 1962;Mason, 1956;Wang, 2019).
To evaluate the accuracy of FSDs of Eq. (6), the wavenumber domain Eq. (4) is treated as a non-linear equation f ðkÞ ¼ 0. The numerically solved wavenumber is compared to the exact wavenumber of FTD presented analytically in Eqs. (7) and (8). Figure 2 demonstrates the accuracy evaluation using an arbitrarily designed pure acoustic velocity c 0 ¼ 2500 m/s, and considering cases of weak, median, and strong attenuation with b ¼(0.010, 0.183, 0.351). The accuracy of the phase velocity is high in general, as the root-mean-square (RMS) error in three cases together is 4.169 m/s. The accuracy of the attenuation depends on the b value. The RMS errors in the attenuation coefficient are 0:0501 Â 10 À3 m À1 for b ¼ 0:010 and 0:3405 Â 10 À3 m À1 for b ¼ 0:190. However, the RMS error in the attenuation coefficient is 6:769 Â 10 À3 m À1 for the extreme case with b ¼ 0:351. This error existed in the most attenuating case is caused by the approximation x=c 0 % k (Eq. 3) used in the transformation from FTD to FSDs.
To investigate the reliability of the FSDs, we compare the analytical solutions of two types of wave equations presented as FTD and FSDs, respectively. Considering a homogeneous case with constant b, the one-dimensional analytical solution for FTD in frequency domain is derived aŝ with whereûðx; xÞ is the wavefield in the frequency domain, andf ðxÞ is the source signature in the frequency domain. For FSDs, the corresponding analytical solution is formed using the Green's function Gðt; k; sÞ as whereûðt; kÞ is the wavefield in the wavenumber domain. The Green's function is where Hðt À sÞ is the Heaviside step function, and An inverse Fourier transform ofûðt; kÞ with respect to the wavenumber k produces the time-space domain wavefield uðt; xÞ. Note that the Green's function in Eq. (12) is presented in terms of a sinc function. We set a model with the velocity of 2500 m/s, and assume the source signature be a Ricker wavelet with the peak frequency of 20 Hz (Wang, 2015). Figure 3 displays the waveform at travel distance 200 m, and demonstrates that two equations match well in general, and only a minor discrepancy exists in the strong attenuation case. The RMS differences corresponding to the cases with b = (0.010, 0.144, 0.190, 0.351) are (5:276 Â 10 À2 , 6:538 Â 10 À2 , 7:295 Â 10 À2 , 22:878 Â 10 À2 ), respectively. This observation is consistent with that shown in Fig. 2.
Noted that the above equation with FSDs (Eq. 6) is derived based on homogeneous media where the viscoelastic parameter b is constant. When forming the FSDs for a general viscoelastic case, bðxÞ is a spatial function. Based on the small perturbation assumption, FSDs is still approximated valid for a general viscoelastic case (Xing & Zhu, 2019;Zhu & Harris, 2014).

The Spatial Filter for the Heterogeneous Media
For the numerical calculation of Eq. (6), we use a pseudo-spectral method to solve the FSD. In practice, the viscoelastic parameter bðxÞ in the heterogeneous media is assumed to be smoothly varied and then the average parameter b is adopted for calculating the FSD: where F x is the Fourier transform with respect to vector x, and F À1 x is the inverse of F x . The wavefield in the space domain and in the wavenumber domain are listed in pairs as the following: uðxÞ $ F x ½uðxÞ; Àr 2 uðxÞ $ k 2 F x ½uðxÞ; The pseudo-spectral method has been used widely in the wave simulation (Carcione, 2010). Therefore, the numerical advantage of the FSD is to overcome the In order to correct the errors caused by the averaging scheme of the FFT implementation, we introduce a correction function f ðbðxÞÞ as a spatial filter to correct the wavenumber as Multiplying either k or k 2 to both sides, we have k 1þbðxÞ ¼ f ðbðxÞÞk 1þb and k 2þbðxÞ ¼ f ðbðxÞÞk 2þb . Therefore, a corrected wave equation which corresponds to Eq. (6) is For constructing the spatial filter, we assume the case of weak attenuation with k Im =jkj ( 1 and make approximationk % k Re . The real wavenumber k Re in Eq. (7) can be expanded to the first order as Then, the spatial filter f ðbðxÞÞ is evaluated at each grid by where x m is the mean frequency. The central-frequency approximation in the last line is made based on an assumption bðxÞ À b ( 1, so that the spatial filter f ðbðxÞÞ is frequency independent and avoids any extra Fourier transform. Thus, this spatial filter does not affect the efficiency of the algorithm, but greatly improve the accuracy in heterogeneous media.

Numerical Examples
In this paper, we present two numerical examples. The first example is for validating the effectiveness of the spatial filter.
We set the constant velocity to be 2500 m/s, and consider a model with the viscoelastic parameter values b ¼ (0.351, 0.190, 0.131). We use the same source signature of 20-Hz Ricker wavelet (Wang, 2015a, b), and calculate the waveforms at distances of 500 m and 2000 m. These two accurate waveforms are plotted in a single trace in Fig. 4 (black solid  curves).
To mimic the approximation in the viscoelastic wave equation, we assume a ''heterogeneous'' model with the average value of b ¼(0.237, 0.152, 0.112). The approximated solutions (in red dots) are overlaid with the accurate trace, as shown in Fig. 4a. Three cases have RMS differences of 23.95 Â10 À2 , 9.09 Â10 À2 , and 4.74 Â10 À2 . These differences are caused mainly due to the phase discrepancy but is also due to the amplitude difference at large b value case.
Adopting the correction with the spatial filter, calculated waveforms (in blue dots) are close to the accurate waveforms, as shown in Fig. 4b. The RMS differences are (2.85 Â10 À2 , 1.05 Â10 À2 , 0.52 Â10 À2 ) for the three cases, respectively. Both the phase and the amplitude are corrected remarkably. This example shows that the correction function can improve the accuracy of waveform simulation in the heterogeneous media, especially in a high-attenuation area.
Next, we validate our scheme in a two-layer model. The model is shown in Fig. 5a. A reference is set by directly solving Eq. (1) using Grünwald-Letnikov expansion (Podlubny, 1999) as follows: where Dt is the time interval. A 20-Hz Ricker wavelet is emitted the centre of the model. The model is discretized into 801 9 801 grids with grid spacings of 2.5 m. This fine spacing is equivalent to 16 nodes per wavelength (k min ¼ v min =ð2:5f p Þ = 40 m). The time step is set as Dt = 0.25 ms, which is also finer than the numerical requirement Dt\Dx=ð ffiffi ffi 2 p v max Þ ¼ 0:59.
Setting fine intervals in spatial and temporal axes is for minimizing discrepancy between Grünwald-Letnikov expansion and the pseudo-spectral method. Figure 5 also displays the wavefield snapshots at 0.35 s of the layered model. The result without correction by the spatial filter (Fig. 5b) shows significant discrepancy from the reference, which proves that the conventional averaging scheme causes errors. However, after correction (Fig. 5c), the accuracy is significantly improved. This example further demonstrates the importance of spatial filter for seismic sim ulation in heterogeneous media. Any remaining weak residual in Fig. 5c is attributed to the transformation from FTD to FSDs.
In the second example, we apply FSDs of Eq. (17) to simulate the wavefield of the Marmousi model. Figure 6 displays the acoustic velocity of the Marmousi model, and the b model. The b model is built based on an analysis of the attenuation versus velocity from a field 3D seismic data in Tarim basin. The model is discretized into 751 9 301 grid points with regular vertical and horizontal grid spacings of 10 m. The source signature is a 20-Hz Ricker wavelet and is emitted at (3800, 150) m. The receivers are located at a depth of 150 m and in a spatial range of 0-7500 m with 10 m spacing. The time step is 1 ms. Figure 7 shows snapshots of the wavefield without attenuation and the wavefield with attenuation at 0.6 s, 0.8 s and 1.0 s, respectively. There are no freesurface reflections from the top boundary and other three boundaries, since a convolutional perfectly matched layer (CPML) absorbing condition is adopted in numerical calculation. However, there are clear reflections from the interfaces within the model, for both cases; this observation indicates the high quality of the simulation with ignorable numerical dispersion. More importantly, these snapshots demonstrate that the attenuating wavefields have a clearly delayed wavefront and reduced amplitude, if compared to its non-attenuation counterparts.
Comparison between the non-attenuating and attenuating shot gathers ( Fig. 8a and b) demonstrates the accumulative effect of attenuation. Moreover, comparison between residuals of the conventional averaging scheme and the corrected scheme with the proposed spatial filter ( Fig. 9a and b) demonstrates the significance of the spatial filter. The residual shown in Fig. 9 is the discrepancy from the Grünwald-Letnikov expansion. Whereas the proposed wave equation is applicable to complex geological model, generated through an empirical formula Vol. 179, (2022) Spatial Filter for the Pseudo-spectral Implementation 2837 models, the reflection events from the interior interfaces are very weak in amplitude. The accurate wave simulation will lead to correct subsurface image from seismic migration, as it leads to correct compensation to the viscoelasticity of the subsurface media.

Conclusions
When transferring the wave equation in the FTD form to the FSDs form, the wave simulation can be implemented in the wavenumber domain. In this paper, we have proposed to insert a frequency-independent correction function into the wave equation as a spatial filter to correct the error caused by the heterogeneity of the model. This spatial filter can be easily implemented and an equation including the correction improves the accuracy of the simulation. Numerical examples have demonstrated that this strategy may properly represent the dissipation effect of the viscoelasticity on waveforms and improve accuracy and maintains the high efficiency of FFT implementation.