Abstract

The research of viscoelastic media is currently a hot topic in the interpretation and processing of seismic data. To accurately simulate the propagation of seismic waves in viscoelastic media, the fractional viscoelastic equation has emerged as an indispensable method. However, solving this equation numerically has proven to be challenging due to the complexity introduced by its fractional Laplacian operators. Recently, deep learning, especially Fourier neural operators (FNO), has shown excellent performance in learning to fast solve partial differential equations. Traditional FNO methods may face crosstalk problems and this make it difficult to achieve satisfactory accuracy when solving the multicomponent fractional order viscoelastic equation. To solve this problem, we introduce a novel approach based on U-net Fourier neural operator (U-FNO). As an enhanced learning method to the traditional FNO-based method, the U-FNO-based method integrates a U-Fourier layer following the standard Fourier layer as a form of regularization, thereby achieving superior prediction accuracy for multicomponent equations. Specifically, both the Fourier layers and U-Fourier layers in U-FNO are trained with the solutions of the equation from previous time steps as inputs. This training process enables the U-FNO to efficiently produce more accurate solutions for subsequent wavefield. Numerical simulations reveal that the U-FNO-based method efficiently learns to solve the fractional viscoelastic wave equation independent of fractional Laplacian operators. Additionally, U-FNO-based method offers superior prediction accuracy in comparison with the traditional FNO-based method.

1. Introduction

Seismic wave forward modelling is essential in the process of seismic data processing (Carcione 1993; Carcione et al.  2002; Liu & Sen 2009a,b; Ren & Liu 2013; Xu et al.  2023). During their propagation through media, seismic waves are affected by attenuation and elasticity: inherent properties of the earth (Robertsson et al.  1994; Alkhalifah 2000; Duveneck & Bakker 2011). The viscoelastic equation, widely recognized for describing the characteristics of viscoelastic media, has garnered considerable attention. Numerous scholars have contributed to the advancement of the fractional viscoelastic wave equation, utilizing constant-Q models to simulate these properties during the propagation of the seismic wavefield (Kalyani et al.  2014). This equation can describe the attenuating effects in an attenuating medium (Li et al.  2016; Zhu 2017; Zhu et al.  2019). However, the solution to the complex equation needs several Fourier transforms because of the fractional Laplacian operators, which are associated with the quality factor Q. Considering the constraints of computing hardware, the necessity for multiple Fourier transforms constitutes a major problem (Yao et al.  2017). Furthermore, to accurately represent the complex propagation of the wavefield, the solution to the viscoelastic equation comprises multiple components, adding to the challenge of solving the equation. Recently, a range of techniques including low-rank decomposition (Sun et al.  2015; Chen et al.  2019), Taylor series expansion (Guo et al.  2016; Zhang et al.  2020), and independent fractional operators (Patnaik et al.  2020; Wang et al.  2022) are introduced into improve the simulation of viscoelastic media. To obtain seismic wavefields of higher frequency and resolution, these methods necessitate finer grid discretization, leading to increased computational costs in forward modelling (Liu & Sen 2011). This limitation often presents a significant challenge in applications such as seismic migration and inversion. Selecting the appropriate solution methods to enhance both accuracy and computational efficiency has emerged as a primary research focus in recent years (Xing & Zhu 2019; Xiong & Guo 2022; Zhou et al.  2023).

The successful implementation of deep learning, particularly through data-driven methods, has inspired numerous scholars to explore suitable approaches for solving partial differential equations (PDEs). A data-driven method for solving PDEs involves learning the mathematical and physical behaviour of these equations from data in a supervised learning framework without relying on complex traditional equations (LeCun et al.  2015; Qu et al.  2022). Recent progress in data-driven-based techniques has dramatically improved computational efficiency, often achieving speeds several orders of magnitude faster than those of traditional solvers. In the realm of data-driven strategies for seismic simulation, two popular methods are widely used. The initial method is the spatial mapping technique, which frequently utilizes convolutional neural networks (CNNs) to apply finite-dimensional operators within its category. These methods have demonstrated remarkable success in quickly and accurately predicting outputs for high-dimensional PDEs (Jiang et al.  2021; Zhong et al.  2023). However, the lack of physical information results in reduced prediction accuracy for this technique (Wang et al.  2017; Smith et al.  2020; Bhattacharya et al.  2021; Wen et al.  2021). The second technique employs the physical information constraint approach, integrating these constraints into a complex network to achieve solutions with high accuracy (Raissi et al.  2019; Zhang et al.  2023). During training, this method (such as physical information neural network (PINN)) automatically satisfies specific physical information using automatic differentiation, leading to enhanced accuracy and improved generalization. Unlike data-driven CNN-based methods, the PINN-based approach falls into the category of wave-equation methods and shares strong resemblances to classical numerical solvers (Moseley et al.  2020; Alkhalifah et al.  2021; Karniadakis et al.  2021; Song et al.  2021; Chen & Ge 2023). The PINN-based method is more suitable for equations that are amenable to solution by finite-difference methods. However, solving this fractional equation with the finite-difference method poses challenges, rendering the PINN-based method unsuitable for equation solving (Wang et al.  2018).

Currently, the Fourier neural operator (FNO) is introduced to solve the problem of PDEs in the supervised learning framework (Li et al.  2020a,b). Differing from the methods previously mentioned, FNO adopts an inductive bias strategy by integrating physical information within the structure of the convolutional blocks (Lu et al.  2019; Grady et al.  2023). The FNO-based method uses a framework that accurately captures the configuration of the wave propagator in the Fourier domain. It is achieved through the use of two fully connected layers, specifically designed for the spatial and wavenumber domains (Kosloff & Baysal 1982; Stoffa et al. 1990). This strategy estimates the mathematical–physical properties of PDEs by establishing the mapping with solutions of the equations at different times and spatial locations (Rashid et al.  2022). Subsequently, the FNO-based method has been widely applied to learning solutions for seismic wave equations (Song & Wang 2022; Zhang et al.  2023). The commonly used training strategy in the time domain is to use the wavefield snapshot at this time step to predict the wavefield snapshot at next time step and establish an appropriate loss function to update the network (Wei & Fu 2022). Recently, by introducing model parameters to modify the FNO network structure, the FNO method has been successfully applied to the rapid inversion of model parameters (Yang et al.  2023; Yin et al.  2023). The FNO-based method has shown great potential for development in the fields of forward modelling and inversion in geophysics. However, its prediction accuracy for complex equations falls short because of the intrinsic regularization effect present in the FNO framework. Meanwhile, solving the viscoelastic equation, especially with multiple component solutions, poses a significant challenge. Recently, U-net has become increasingly popular as a deep learning architecture for data analysis, especially in image processing (Long et al.  2015). It uses convolutional layers to extract features and has succeeded in tasks like image classification and object detection (Krizhevsky et al.  2012; He et al.  2016). Combining the benefits of both FNO and U-net, researchers have proposed the U-net Fourier neural operator (U-FNO) and solved the problems with complex equations, such as multiphase flow equations (Wen et al.  2022). Their research declares that the U-FNO-based method outperforms the traditional FNO in improving the accuracy of complex PDE solutions. Currently, U-FNO has yet to be applied in the field of geophysics.

We examine the performance of FNOs in learning viscoelastic wave equations and compare the predictive performance between the standard FNO-based method and the U-FNO-based method. The U-FNO method combines the advantages of both FNO- and U-net-based methods. The U-net, employed as a regularization term, effectively addresses the crosstalk problem in traditional FNO methods for multicomponent equations and achieves high-accuracy solutions for multicomponent fractional viscoelastic equations. This provides a foundation for achieving high-precision inversion by modifying the network.

This paper first introduces the fractional viscoelastic equation. Subsequently, we offer an overview of the FNO framework used for seismic modelling. Next, the theory of the U-FNO architecture for seismic modelling is introduced. Finally, we assess the precision and performance of our proposed method by a homogeneous model and a partial Hess model.

2. Methodology

In this section, we detail the fractional viscoelastic equation along with its solving techniques, discuss the structure and problem of the FNO, and introduce the U-FNO framework as an innovative solution to these problems.

2.1. The introduction of fractional viscoelastic equation

The 2D fractional viscoelastic equation, derived from the constant-Q model, is often formulated as follows (Zhu & Carcione 2014):

(1)

where v=(vx,vz) and t=(τxx,τzx,τzz) represent the particle velocity and stress components, respectively, ρ denotes the density, and Dξ can be described as follows:

(2)

where ξ means the parameters related to the P- and S-waves individually; γξ=arctan(1/Qξ)/π represents the fractional operators, which are associated with the quality factor Qξ; 2 denotes the Laplacian operator, aξ and bξ are the intermediate variables that can be described as follows:

(3)

where ω0=2πf0 stands for the reference angular frequency and Vξ=(Vp0,Vs0) represent the reference P-velocity and reference S-velocity.

Owing to the involvement of two fractional Laplacian operators (2)γξ and (2)γξ0.5, which are dependent on the quality factor Q and exhibit variations across different spatial positions, solving them using the conventional finite-difference method remains a substantial challenge. Currently, several methods have been proposed to achieve efficient and accurate solutions for fractional Laplacian operators, enabling rapid and precise propagation of viscoelastic wavefields. In this paper, the staggered grid pseudo-spectral (SGPS) method is employed as a widely used technique for solving the viscoelastic equation to obtain the training and testing data. To enable the SGPS method to solve variable fractional equations of complex models, we approximate the variable Q value with an average Q value. In solving Equation (1), the solution frequently manifests as particle velocity, represented by v=(vx,vz). To obtain this solution, we further rewrite Equation (1) as follows:

(4)

where t and t1 represent this time step and the previous time step, respectively. This can be further simplified to the following form:

(5)

where G is a function representing Equation (4). From Equation (5), to predict the wavefield snapshot at this time step, we only need the wavefield at previous time step. To simplify the subsequent use of symbols we let vt and vt1 stand for v(t) and v(t1), respectively.

2.2. The review of FNO architecture for seismic modelling

Currently, an innovative method known as the FNO is emerged, which entails the direct parameterization of the integral kernel within Fourier domain. Its main objective is to build a mapping in an infinite-dimensional space with a limited set of an input–output dataset, a critical task in geophysics for modelling complex wave propagation. This can be accomplished by employing convolution with functions related to low wavenumbers to capture global features, and then applying an activation function to these features to transition back to the high wavenumber model. This process is crucial for modelling seismic wave propagation. The exceptional capability of the FNO allows it to approximate functions characterized by Fourier modes. Its architecture is comprised of multiple Fourier layers, with each layer hosting two fully connected layers in both the spatial and wavenumber domains, as shown in Fig. 1b. The fully connected layers in the spatial domain, which includes a local transform W, plays a pivotal role as the trainable phase-screen compensator, enabling a model to accommodate local variations. Meanwhile, the fully connected layer within the wavenumber domain integrates a sequence of transformations: a Fourier transform F, followed by a linear transformation R, and concludes with an inverse Fourier transform F1; It serves as a non-local spatial convolution operator and functions as a wavenumber filter within the process of phase-shift. We introduce a four-layer FNO architecture designed for learning to solve the fractional viscoelastic equation, illustrated in Fig. 1a. The FNO layer receives two-component wavefields at a partial time step as input and predicts the two-component wavefields for the subsequent time step. The mathematical representation of this learning process is as follows:

(6)

where i1 represents the previous time step and FC1 is a full connection layer that elevates the data to a higher-dimensional channel. Equation (6) reveals the data flow between the FNO time-space and time-wavenumber domains. The predicted wavefield at this time step i is:

(7)

where FC2 includes the two full connection layers and a ReLU function, and it make data back to the target dimension. Equations (6) and (7) can also be described as:

(8)

where O stands for a nonlinear map operator.

Training and prediction processes by FNO (a) and U-FNO (c) and the architectures of Fourier layers (b) and U-Fourier layers (d).
Figure 1.

Training and prediction processes by FNO (a) and U-FNO (c) and the architectures of Fourier layers (b) and U-Fourier layers (d).

By comparing Equation (5) with Equation (8), it becomes clear that the method based on the FNO seeks to refine the nonlinear mapping operator O to capture the mathematical–physical characteristics of the function G. This objective can be attained with optimizing a loss function that evaluates the difference between the predicted O(vt1) and true wavefields G(vt1) through the L2 norm, a widely accepted metric for quantifying the difference between two datasets. This process allows the FNO to capture the complex interactions and dynamics within the seismic wavefield, thereby enhancing its predictive power and accuracy in solving the fractional viscoelastic equation. Unfortunately, when the FNO-based method solves multicomponent equations, crosstalk occurs between the two components, affecting the network's accuracy in solving equations.

2.3. The theory of U-FNO architecture for seismic modelling

U-FNO is introduced to enhance the optimization of FNO. It consists of the U-Fourier layer (Fig. 1d) and the Fourier layer (Fig. 1b). By incorporating a U-net architecture, the U-Fourier layer improves its capability for parameterizing high-frequency information in the wavefield through local convolutional kernels, distinguishing it from the traditional Fourier layer. This addition increases the flexibility and robustness of FNO, making it more effective in solving the fractional viscoelastic wave equation and in handling complex seismic wavefield modelling. The training process in wavefield modelling can be revised by employing the U-FNO-based method, as shown in Fig. 1c. Specifically, the two-component wavefield with partial time steps is sequentially processed through the FNO layer followed by the U-FNO layer. The process of U-FNO, which is aimed at obtaining the solution of the fractional viscoelastic equation, is as follows:

  • The U-FNO employs the wavefield snapshots from the previous several time steps (Nt) as input, which is then elevated to a higher-dimensional channel domain through a fully connected layer (FC1), leading to a promoted snapshots of wavefield (V).

  • This promoted snapshots of wavefield (V) is further refined with some Fourier layers. In each Fourier layer, the promoted wavefield is transformed into the wavenumber domain by a sequence of operations: Fourier transform F, linear transform R, and inverse Fourier transform F1. Meanwhile, this promoted wavefield is linearly transformed by local transform (W) in the spatial domain. Moreover, the outputs of these two parts are combined by activation function ReLU to obtain the output from this Fourier layer. Subsequently, this output is input into the next Fourier layer.

  • The final output from these Fourier layer (FoutL) serves as the input for several U-Fourier layers. In each U-Fourier layer, the input is processed to obtain (UFout) not only through two components of the conventional FNO but also via the U-net (UNET).

  • The final output, which is obtained by the last U-Fourier layer (UFoutN), is then mapped back to the target dimension using a fully connected layer (FC2). This network includes two fully connected layers and a ReLU function, which enables it to predict the wavefield at the current time.

Similar to Equation (6), the learning processes of FNO and U-FNO can be mathematically represented as follows:

(9)

The final predicted wavefield at this time step (i) is.

(10)

Through training U-FNO with snapshots of solutions to this complex equation, the resulting U-FNO becomes a PDE solver capable of generating accurate solutions for Equation (1). By comparing Equations (9) and (6), it is demonstrated that the features of the wavefield snapshots extracted by the U-net in the U-FNO method can be used as regularization terms to mitigate the crosstalk problem between the two component datasets. The layer number of Fourier (L) and U-Fourier (N) is often tuned based on the characteristics and requirements of the specific problem being solved. After many experiments, we find that employing two Fourier layers followed by two U-Fourier layers can obtain the optimal performance for solving a fractional viscoelastic equation. The input and training approach utilized in this method play a critical role in the training process. Figure 2 offers a schematic depiction of the input and output within the U-FNO (taking snapshots Nt = 3 with as an example). During training, we utilize Nt true snapshots in a recursive method to generate the subsequent one snapshot. Meanwhile, in the prediction process, we employ Nt snapshots, including snapshots of the true and predicted wavefields, to generate the subsequent wavefield snapshot. It is worth noting that in the recursive prediction process, the true wavefield snapshots are decreasing.

Schematic description of input and output in training and prediction (number of snapshots is three).
Figure 2.

Schematic description of input and output in training and prediction (number of snapshots is three).

3. Numerical examples

In this section, the capability of U-FNO to predict seismic wavefield propagation in viscoelastic media is demonstrated by a homogeneous model and a partial Hess model. The SGPS-based method, in conjunction with PML boundary conditions, is employed to solve the fractional viscoelastic wave equation for making dataset. These methodologies allow us to generate datasets for training, validation, and prediction. Importantly, informed by insights from previous research and extensive experiments, we set the dimensions of the Fourier and U-Fourier layers to 12 × 20. Additionally, we use structural similarity (SSIM) and signal-to-noise ratio (SNR) (see Appendix  A) to evaluate the quality of the predicted wavefield snapshots derived from both the homogeneous and partial Hess models.

3.1. The test on a homogeneous model

For the first test, a 2D homogeneous model is employed. This model is designed with a grid size of 100 × 100 and a grid spacing of 10 m. In Table 1, the parameters for velocity and the quality factor Q associated with P- and S-waves are detailed. We use a Ricker wavelet at 35 Hz with an interval of 1 ms, and a maximum sampling time of 0.24 s as the source. Figure 3 provides a visual representation of the distribution of training, validation, and prediction sets. This figure depicts the even distribution of 625 sources across the velocity model, each spaced at 40 m. From shots 100 to 500 (highlighted in yellow), the 314th shot, located at the model's centre and marked by a red star, is selected to evaluate the prediction performance. Here, 80% of the yellow region is randomly allocated to the training dataset, while the remaining 20% is allocated to the validation dataset. These datasets are compiled with a maximum time interval limited to 0.18 s, and particle velocity snapshots are resampled every 10 ms.

A visual depiction of the distribution of the training, validation, and prediction sets.
Figure 3.

A visual depiction of the distribution of the training, validation, and prediction sets.

Table 1.

Model parameters of the 2D homogeneous model.

SuperscriptVp(m/s)Vs(m/s)QpQs
Parameter350025004030
SuperscriptVp(m/s)Vs(m/s)QpQs
Parameter350025004030
Table 1.

Model parameters of the 2D homogeneous model.

SuperscriptVp(m/s)Vs(m/s)QpQs
Parameter350025004030
SuperscriptVp(m/s)Vs(m/s)QpQs
Parameter350025004030

Respectively training by FNO-based and U-FNO-based approaches, we present their predictions for the vx and vz components in Figs. 4 and 5. Figures 4b, d, f, and h and 5b, d, f, and h demonstrate the predictive capabilities of both FNO-based and U-FNO-based approaches in simulating the reflected wavefield, direct wavefield, and the boundaries condition. Our research evaluates the impact of the number with wavefield snapshots on prediction. Through a comparison of Fig. 4b and f with 4d and h, a distinct improvement in P-wave amplitude accuracy emerges with the increment in the number of particle velocity inputs within both the FNO and U-FNO. Meanwhile, the difference (Fig. 4c, e, g, and i) between true wavefield snapshots and the corresponding predicted wavefield snapshots provides additional confirmation of this observation. During the subsequent analysis of Fig. 4b and d and 4f and h, employing an equal number of inputting wavefield snapshots, it becomes evident that the U-FNO amplifies prediction accuracy in comparison to the FNO method. Furthermore, it is worth noting that utilizing a separate wavefield snapshot as the input for the FNO method leads to inadequate prediction accuracy of P-wavefield snapshots. By contrast, U-FNO-based methods can achieve high-precision predictions of wavefields using just a single wavefield snapshot as the input (Fig. 4f). Figure 5 shows the prediction of particle velocity vz, similar with the earlier observations made regarding particle velocity. Simultaneously, through a comparison of the difference in Figs. 4 and 5, it becomes apparent that the predictive accuracy of vz is relatively lower when contrasted with vx, especially in the FNO-based method. This is because the vz component of particle velocity experiences more frequent event polarity reversals compared to the vx component, resulting in relatively poorer event continuity. Consequently, the vz component is more challenging to predict than the vx component. However, using the U-FNO-based method can further improve the prediction accuracy of the vz component.

Snapshots of particle velocity ${{v}_x}$ in the homogeneous model: true particle velocity (a); predicted particle velocity by FNO with inputting one-particle velocity (b) and three-particle velocity (d); predicted particle velocity by U-FNO with inputting one-particle velocity (c) and three-particle velocity (e); the difference between (b) and (a), (c) and (a), (d) and (a), (e) and (a), respectively (f–i).
Figure 4.

Snapshots of particle velocity vx in the homogeneous model: true particle velocity (a); predicted particle velocity by FNO with inputting one-particle velocity (b) and three-particle velocity (d); predicted particle velocity by U-FNO with inputting one-particle velocity (c) and three-particle velocity (e); the difference between (b) and (a), (c) and (a), (d) and (a), (e) and (a), respectively (f–i).

Snapshots of particle velocity ${{v}_z}$ in the homogeneous model: true particle velocity (a); predicted particle velocity by FNO with inputting one-particle velocity (b) and three-particle velocity (d); predicted particle velocity by U-FNO with inputting one-particle velocity (c) and three-particle velocity (e); the difference between (b) and (a), (c) and (a), (d) and (a), (e) and (a), respectively (f–i).
Figure 5.

Snapshots of particle velocity vz in the homogeneous model: true particle velocity (a); predicted particle velocity by FNO with inputting one-particle velocity (b) and three-particle velocity (d); predicted particle velocity by U-FNO with inputting one-particle velocity (c) and three-particle velocity (e); the difference between (b) and (a), (c) and (a), (d) and (a), (e) and (a), respectively (f–i).

To provide further validation for the effectiveness of our method, we enhance our analysis by vertical profiles of particle velocity snapshots. Figure 6 illustrates the particle velocity vx at CDP of 0.5 km, and Fig. 7 showcases the particle velocity vz at CDP of 0.3 km. After making a comprehensive comparison between Figs. 6 and 7, it becomes clear that both the U-FNO-based method and the number of inputting particle velocities enhance the precision in predicting particle velocities (vx, vz). Remarkably, even with a single input particle velocity, the predictions generated using the U-FNO exhibit a better similarity with the true wavefield. This shows the exceptional capacity of the U-FNO to provide accurate and reliable predictions.

Comparison of the vertical profiles in particle velocity ${{v}_x}$ at 0.5 km: vertical profiles of true particle velocity and predicted particle velocity by FNO with inputting one-particle velocity (a) and inputting three-particle velocity (c); vertical profiles of true particle velocity and predicted particle velocity by U-FNO with inputting one-particle velocity (b) and inputting three-particle velocity (d).
Figure 6.

Comparison of the vertical profiles in particle velocity vx at 0.5 km: vertical profiles of true particle velocity and predicted particle velocity by FNO with inputting one-particle velocity (a) and inputting three-particle velocity (c); vertical profiles of true particle velocity and predicted particle velocity by U-FNO with inputting one-particle velocity (b) and inputting three-particle velocity (d).

Comparison of the vertical profiles in particle velocity ${{v}_z}$ at 0.5 km: vertical profiles of true particle velocity and predicted particle velocity by FNO with inputting one-particle velocity (a) and inputting three-particle velocity (c); vertical profiles of true particle velocity and predicted particle velocity by U-FNO with inputting one-particle velocity (b) and inputting three-particle velocity (d).
Figure 7.

Comparison of the vertical profiles in particle velocity vz at 0.5 km: vertical profiles of true particle velocity and predicted particle velocity by FNO with inputting one-particle velocity (a) and inputting three-particle velocity (c); vertical profiles of true particle velocity and predicted particle velocity by U-FNO with inputting one-particle velocity (b) and inputting three-particle velocity (d).

Furthermore, Table 2 presents the computed SSIM and SNR values, providing additional support for the aforementioned observation. Simultaneously, a clear trend emerges from the information in Table 2 is that the prediction of vz poses a more complex prediction, compared with the prediction of vx within both the FNO-based and U-FNO-based methods. Figure 8 illustrates the contrast in training and validation losses between FNO and U-FNO, utilizing one (Fig. 8a) and three (Fig. 8b) input wavefield snapshots. Compared with the FNO-based method, the U-FNO-based approach reveals better convergence. Moreover, as the number of input wavefields increases, both the U-FNO and FNO methods exhibit substantial convergence enhancements for the uncomplicated homogeneous model.

Training losses and validation losses by inputting one snapshot of particle velocity (a) and three snapshots of particle velocity (b).
Figure 8.

Training losses and validation losses by inputting one snapshot of particle velocity (a) and three snapshots of particle velocity (b).

Table 2.

SSIM and SNR of the prediction by FNO and U-FNO with the homogeneous model.

 FNOU-FNO
SuperscriptOne snapshotThree snapshotsOne snapshotThree snapshots
SSIMvx0.30750.65040.69440.8805
 vz0.29240.61130.63730.8773
SNRvx13.7518.9619.6828.57
 vz9.7116.1916.6225.10
 FNOU-FNO
SuperscriptOne snapshotThree snapshotsOne snapshotThree snapshots
SSIMvx0.30750.65040.69440.8805
 vz0.29240.61130.63730.8773
SNRvx13.7518.9619.6828.57
 vz9.7116.1916.6225.10
Table 2.

SSIM and SNR of the prediction by FNO and U-FNO with the homogeneous model.

 FNOU-FNO
SuperscriptOne snapshotThree snapshotsOne snapshotThree snapshots
SSIMvx0.30750.65040.69440.8805
 vz0.29240.61130.63730.8773
SNRvx13.7518.9619.6828.57
 vz9.7116.1916.6225.10
 FNOU-FNO
SuperscriptOne snapshotThree snapshotsOne snapshotThree snapshots
SSIMvx0.30750.65040.69440.8805
 vz0.29240.61130.63730.8773
SNRvx13.7518.9619.6828.57
 vz9.7116.1916.6225.10

3.2. The test on a partial Hess model

To evaluate the adaptability of U-FNO to complex models, we select a partial Hess model for assessing the accuracy of the predicted particle velocities (vx, vz). Figure 9 displays the P- and S-wave velocities, along with the quality factors Q. This partial Hess model is divided into a 100 × 100 grid with an 8-m grid interval. A Ricker wavelet with main frequency of 35 Hz serves as the source wavelet. Its maximum sampling time and time interval are 0.25 s and 0.5 ms, respectively. To improve prediction accuracy, we set the inputting snapshot to 5. The particle velocity snapshots are taken every 3 ms for resampling. The settings of the training dataset and the validation dataset are similar to those of the homogeneous model (Fig. 3). These datasets have a maximum sampling time of 0.13 s.

Partial Hess model: (a) Vp; (b) Vs; (c) Qp; and (d) Qs.
Figure 9.

Partial Hess model: (a) Vp; (b) Vs; (c) Qp; and (d) Qs.

Figures 10 and 11 provide comparisons of predicted particle velocity snapshots at 0.15 s, a time step that extends beyond the training period respectively. In Figure 10b and d, both methods reconstruct the particle velocity snapshots of the complex model. However, on comparing Fig. 10c with e, it becomes evident that the difference between the predicted particle velocity snapshots vx obtained through the U-FNO and the true particle velocity is reduced compared to the difference observed in the FNO method. This observation indicates that the U-FNO approach demonstrates enhanced accuracy in predicting the particle velocity vx. for complex models. This trend is similarly noticeable in the prediction of the particle velocity vz. On comparing Figs. 10 and 11, it becomes evident that predicting the particle velocity vz is more challenging than predicting vx for complex models. Nevertheless, U-FNO exhibits superior performance in predicting the particle velocity vz. To achieve a clear comparison of the accuracy for particle velocity snapshots within this complex model, we have extracted vertical profiles of particle velocity vz at 0.6 km and particle velocity vx at 0.7 km, as illustrated in Fig. 12. When Fig. 12 parts a–c are compared with Fig. 12 parts b–d, it becomes clear that the particle velocity snapshots (vz, vx) predicted by U-FNO bear a stronger resemblance to the true particle velocity snapshots (vz, vx) in comparison to the predictions made by FNO. We also provide a comparison of the training and validation losses for both FNO and U-FNO using the partial Hess model in Fig. 13. For complex models, U-FNO also demonstrates enhanced generalization and convergence capabilities in comparison to FNO.

Snapshots of particle velocity ${{v}_x}$ in a partial Hess model: true particle velocity (a); predicted particle velocity by FNO-based method (b) and U-FNO-based method (c); the differences between (b) and (a), (c) and (a) (d–e).
Figure 10.

Snapshots of particle velocity vx in a partial Hess model: true particle velocity (a); predicted particle velocity by FNO-based method (b) and U-FNO-based method (c); the differences between (b) and (a), (c) and (a) (d–e).

Snapshots of particle velocity ${{v}_z}$ in partial Hess model: true particle velocity (a); predicted particle velocity by FNO-based method (b) and U-FNO-based method (c); the differences between (b) and (a), (c) and (a) (d–e).
Figure 11.

Snapshots of particle velocity vz in partial Hess model: true particle velocity (a); predicted particle velocity by FNO-based method (b) and U-FNO-based method (c); the differences between (b) and (a), (c) and (a) (d–e).

Comparison of the vertical profiles in true and predicted particle velocity ${{v}_x}$ by FNO (a) and U-FNO (b); comparison of the vertical profiles in true and predicted particle velocity ${{v}_z}$ by FNO (c) and U-FNO (d).
Figure 12.

Comparison of the vertical profiles in true and predicted particle velocity vx by FNO (a) and U-FNO (b); comparison of the vertical profiles in true and predicted particle velocity vz by FNO (c) and U-FNO (d).

Training and validation losses of epochs for partial Hess model.
Figure 13.

Training and validation losses of epochs for partial Hess model.

The differences observed in Figs. 14c and f and 15c and f indicate that as seismic waves propagate, both methods experience errors after a long period of propagation. However, compared to traditional FNO-based methods, our proposed method maintains commendable predictive performance. Simultaneously, we further evaluate the training performance of our proposed method on wavefields with source locations distant from the training set. We select the wavefield snapshots from the 614th source location as the prediction set. As shown in Fig. 16, the wavefield of the 614th source is already distant from the training set, and its prediction performance can well illustrate the generalization ability of our method. From the wavefield snapshot predicted at 0.17 s for the 614th shot (Fig. 17), it can be observed that even when far from the training region, our proposed method demonstrates good prediction performance. This indicates that our method has better generalization. U-FNO-based method exhibits better predictive performance in component (vx,vz) than FNO-based method. From Table 3, SSIM and PSNR demonstrate that as seismic waves propagate over a long period, the prediction accuracy of both FNO and U-FNO decreases to varying degrees. However, the U-FNO-based method achieves better performance in prediction. In addition, the SSIM and SNR values at the 614th frame at 0.17s indicate that the U-FNO-based method has higher accuracy than the FNO-based method when away from the training data area. Observations in Table 3 indicate that our proposed method has better generalization ability and can make better predictions for data that is longer in duration and farther away from the training area.

Snapshots of particle velocity in 314th shot at 0.17 s: true particle velocity ${{v}_z}$(a) and ${{v}_x}$(f), predicted particle velocity ${{v}_z}$(b) and ${{v}_x}$(g) by FNO-based method, the difference in ${{v}_z}$(c) and ${{v}_x}$(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity ${{v}_z}$(d) and ${{v}_x}$(i) by U-FNO-based method, the difference in ${{v}_z}$(e) and ${{v}_x}$(j) between true particle velocity and predicted particle velocity by U-FNO-based method.
Figure 14.

Snapshots of particle velocity in 314th shot at 0.17 s: true particle velocity vz(a) and vx(f), predicted particle velocity vz(b) and vx(g) by FNO-based method, the difference in vz(c) and vx(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity vz(d) and vx(i) by U-FNO-based method, the difference in vz(e) and vx(j) between true particle velocity and predicted particle velocity by U-FNO-based method.

Snapshots of particle velocity in 314th shot at 0.19 s: true particle velocity ${{v}_z}$(a) and ${{v}_x}$(f), predicted particle velocity ${{v}_z}$(b) and ${{v}_x}$(g) by FNO-based method, the difference in ${{v}_z}$(c) and ${{v}_x}$(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity ${{v}_z}$(d) and ${{v}_x}$(i) by U-FNO-based method, the difference in ${{v}_z}$(e) and ${{v}_x}$(j) between true particle velocity and predicted particle velocity by U-FNO-based method.
Figure 15.

Snapshots of particle velocity in 314th shot at 0.19 s: true particle velocity vz(a) and vx(f), predicted particle velocity vz(b) and vx(g) by FNO-based method, the difference in vz(c) and vx(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity vz(d) and vx(i) by U-FNO-based method, the difference in vz(e) and vx(j) between true particle velocity and predicted particle velocity by U-FNO-based method.

A visual depiction of the distribution of 614th shot prediction dataset.
Figure 16.

A visual depiction of the distribution of 614th shot prediction dataset.

Snapshots of particle velocity in 614th shot at 0.17 s: True particle velocity ${{v}_z}$(a) and ${{v}_x}$(f), predicted particle velocity ${{v}_z}$(b) and ${{v}_x}$(g) by FNO-based method, the difference in ${{v}_z}$(c) and ${{v}_x}$(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity ${{v}_z}$(d) and ${{v}_x}$(i) by U-FNO-based method, the difference in ${{v}_z}$(e) and ${{v}_x}$(j) between true particle velocity and predicted particle velocity by U-FNO-based method.
Figure 17.

Snapshots of particle velocity in 614th shot at 0.17 s: True particle velocity vz(a) and vx(f), predicted particle velocity vz(b) and vx(g) by FNO-based method, the difference in vz(c) and vx(h) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity vz(d) and vx(i) by U-FNO-based method, the difference in vz(e) and vx(j) between true particle velocity and predicted particle velocity by U-FNO-based method.

Table 3.

SSIM and SNR of the prediction by FNO and U-FNO with the partial Hess model.

   
 SSIMSNR
   
 vxvzvxvz
   
SuperscriptFNOU-FNOFNOU-FNOFNOU-FNOFNOU-FNO
314th at 0.15 s0.88770.93660.86680.914532.8137.3729.134.06
314th at 0.17 s0.88410.93570.86520.913130.135.8726.5932.96
314th at 0.19 s0.88330.93420.86320.912728.1334.9124.5931.33
614th at 0.17 s0.64630.72480.60630.716326.8529.4723.7327.3
   
 SSIMSNR
   
 vxvzvxvz
   
SuperscriptFNOU-FNOFNOU-FNOFNOU-FNOFNOU-FNO
314th at 0.15 s0.88770.93660.86680.914532.8137.3729.134.06
314th at 0.17 s0.88410.93570.86520.913130.135.8726.5932.96
314th at 0.19 s0.88330.93420.86320.912728.1334.9124.5931.33
614th at 0.17 s0.64630.72480.60630.716326.8529.4723.7327.3
Table 3.

SSIM and SNR of the prediction by FNO and U-FNO with the partial Hess model.

   
 SSIMSNR
   
 vxvzvxvz
   
SuperscriptFNOU-FNOFNOU-FNOFNOU-FNOFNOU-FNO
314th at 0.15 s0.88770.93660.86680.914532.8137.3729.134.06
314th at 0.17 s0.88410.93570.86520.913130.135.8726.5932.96
314th at 0.19 s0.88330.93420.86320.912728.1334.9124.5931.33
614th at 0.17 s0.64630.72480.60630.716326.8529.4723.7327.3
   
 SSIMSNR
   
 vxvzvxvz
   
SuperscriptFNOU-FNOFNOU-FNOFNOU-FNOFNOU-FNO
314th at 0.15 s0.88770.93660.86680.914532.8137.3729.134.06
314th at 0.17 s0.88410.93570.86520.913130.135.8726.5932.96
314th at 0.19 s0.88330.93420.86320.912728.1334.9124.5931.33
614th at 0.17 s0.64630.72480.60630.716326.8529.4723.7327.3

3.3. The test on a partial Marmousi model

To verify the applicability of the proposed method for complex models, we apply the partial Marmousi model in this section. Figure 18 presents the P- and S-wave velocities of the partial Marmousi model, along with the quality factor Q. All these models are divided into a 100 × 100 grid with grid intervals of 10 m. A Ricker wavelet, with a central frequency of 30 Hz, is used as the source wavelet. The maximum sampling duration and the time interval are 0.2 s and 0.5 ms, respectively. In this model test, we load the seismic source simultaneously on both the components (vx, vz) to verify the reason for the low prediction accuracy of the vz component. The input snapshots are set to 5. The particle velocity snapshots are captured every 3 ms for resampling. The settings for the training dataset and validation dataset are shown in Fig. 3. These datasets have a maximum sampling time of 0.13 s.

Partial Marmousi model: (a) Vp; (b) Vs; (c) Qp; and (d) Qs.
Figure 18.

Partial Marmousi model: (a) Vp; (b) Vs; (c) Qp; and (d) Qs.

Figures 19 and 20 compare the predicted particle velocity snapshots at 0.17 s, which is a time step beyond the training period. From Fig. 19b and d, it can be observed that both FNO and U-FNO methods can be effectively applied to complex models. However, by comparing Fig. 19c and e, it is evident that the differences obtained using the U-FNO method are significantly smaller than those using the FNO method. This indicates that our proposed method has higher prediction accuracy. Similar observations can also be seen in Fig. 20. It is worth noting that the difference in the vz component of this model does not significantly increase compared to the vx component. This is because when the seismic source loads both the components (vx,vz) simultaneously, the continuity of events in the vz component is similar to that in the vx component, and there is no significant polarity reversal. In addition, we also use the SSIM and SNR in Table 4 to verify the superiority of our proposed method. As shown in Table 4, the SSIM and SNR of our proposed method for multicomponent prediction are both higher than those of the traditional FNO method. This indicates that our proposed method has higher prediction accuracy.

Snapshots of particle velocity in 314th shot at 0.17 s: True particle velocity ${{v}_x}$(a), predicted particle velocity ${{v}_x}$(b) by FNO-based method, the difference in ${{v}_x}$(c) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity ${{v}_x}$(d) by U-FNO-based method, the difference in ${{v}_x}$(e) between true particle velocity and predicted particle velocity by U-FNO-based method.
Figure 19.

Snapshots of particle velocity in 314th shot at 0.17 s: True particle velocity vx(a), predicted particle velocity vx(b) by FNO-based method, the difference in vx(c) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity vx(d) by U-FNO-based method, the difference in vx(e) between true particle velocity and predicted particle velocity by U-FNO-based method.

Snapshots of particle velocity in 314th shot at 0.17 s: true particle velocity ${{v}_z}$(a), predicted particle velocity ${{v}_z}$(b) by FNO-based method, the difference in ${{v}_z}$(c) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity ${{v}_z}$(d) by U-FNO-based method, the difference in ${{v}_z}$(e) between true particle velocity and predicted particle velocity by U-FNO-based method.
Figure 20.

Snapshots of particle velocity in 314th shot at 0.17 s: true particle velocity vz(a), predicted particle velocity vz(b) by FNO-based method, the difference in vz(c) between true particle velocity and predicted particle velocity by FNO-based method, predicted particle velocity vz(d) by U-FNO-based method, the difference in vz(e) between true particle velocity and predicted particle velocity by U-FNO-based method.

Table 4.

SSIM and SNR of the prediction by FNO and U-FNO with the partial Marmousi model.

SuperscriptFNOU-FNO
SSIMvz0.84900.8994
 vx0.84630.9012
SNRvz23.5129.67
 vx24.9930.42
SuperscriptFNOU-FNO
SSIMvz0.84900.8994
 vx0.84630.9012
SNRvz23.5129.67
 vx24.9930.42
Table 4.

SSIM and SNR of the prediction by FNO and U-FNO with the partial Marmousi model.

SuperscriptFNOU-FNO
SSIMvz0.84900.8994
 vx0.84630.9012
SNRvz23.5129.67
 vx24.9930.42
SuperscriptFNOU-FNO
SSIMvz0.84900.8994
 vx0.84630.9012
SNRvz23.5129.67
 vx24.9930.42

3.4. Computational performance

For a more comprehensive comparison of computational performance, we have detailed the computational costs in Table 5, with all training, testing, and prediction processes conducted on a GPU (GeForce RTX 3060). The computational performance of the traditional method is not included in Table 4 due to its implementation on the GPU using C code, rendering comparisons with CUDA-based Python implementations less meaningful without adequate context. Notably, the computational efficiency of the FNO method in generating wavefields has shown significant acceleration compared to traditional methods once the network is adeptly trained. This efficiency improvement is quantified as an enhancement of two to three orders of magnitude in computational speed. In our key comparison, we evaluated the computational performance of dataset prediction on the GPU for both FNO and U-FNO methods. We found that the computational cost for U-FNO is 0.61 times higher than that of the FNO method. This increase can be ascribed to the integration of an additional U-net layer in conjunction with the conventional Fourier layers. Nonetheless, U-FNO demonstrates a notable improvement in computational performance when contrasted with the traditional SGPS method.

Table 5.

Computational performance of FNO and U-FNO.

MethodNumber of trainable parametersTime (s)
FNO926 4170.00 353
U-FNO1 156 3370.00 571
MethodNumber of trainable parametersTime (s)
FNO926 4170.00 353
U-FNO1 156 3370.00 571
Table 5.

Computational performance of FNO and U-FNO.

MethodNumber of trainable parametersTime (s)
FNO926 4170.00 353
U-FNO1 156 3370.00 571
MethodNumber of trainable parametersTime (s)
FNO926 4170.00 353
U-FNO1 156 3370.00 571

4. Conclusions

Our research aims to develop a more accurate surrogate model to approximate the solution of viscoelastic wave equations while maintaining an acceptable computational cost. Compared with the traditional FNO-based method, our method introduces U-Fourier layers subsequent to the standard Fourier layers. This improvement effectively approximates the nonlinear mapping operator that governs solution of the equation through data-driven training, thereby markedly improving the accuracy solution of equations. Notably, our approach retains the benefits of the conventional FNO-based method, including obviating the need for the fractional Laplace operator and enhancing computational efficiency. Numerical examples show that our method exhibits better performance in approximating solutions to fractional viscoelastic equations compared to traditional FNO-based methods. Additionally, it can prove significant robustness and generalization ability by predicting solutions for multicomponent fractional viscoelastic equations, even beyond the confines of the dataset. The proposed method demonstrates significant advantages in approximating wave-equation solutions, showing potential for applications that require multiple iterations of simulation. In the future, by optimizing the network and integrating the velocity parameter model, this method could be further enhanced and more effectively applied to the field of parameter inversion.

Acknowledgements

Many colleagues have helped with suggestions for improving this papers.

Conflict of interest statement

The authors declared that they do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Funding

This work is supported by the National Natural Science Foundation of China (grant nos. 42274147 and 41874144).

Data availability

Data associated with this research are available and can be obtained by contacting the corresponding author.

Appendix A

SSIM function is a quantitative metric used to evaluate the similarity between the true snapshot (x) and the predicted snapshot (y). This evaluation is based on factors such as the local mean values (μx and μy), standard deviations (σx2 and σy2), cross-covariances (σxy) of the two models. c1 and c2 are regularization constants, with values of c1=6.5025×108and c2=5.85225×107, respectively. Meanwhile, SNR between the true snapshot (x) and the predicted snapshot (y) quantifies the ratio of signal power to noise power within the models. They are described as.

(A-1)

and

(A-2)

Higher SSIM and SNR values indicate greater similarity between the true and predicted wavefield snapshots.

References

Alkhalifah
 
T
.,
2000
.
An acoustic wave equation for anisotropic media
.
Geophysics
,
65
:
1239
50
.

Alkhalifah
,
T
,
Song
 
C
,
Waheed
 
UB
,
Hao
 
Q
.,
2021
.
Wavefield solutions from machine learned functions constrained by the Helmholtz equation
.
Artif Intell Geosci
,
2
:
11
9
.

Bhattacharya
 
K
,
Hosseini
 
B
,
Kovachki
 
NB
,
Stuart
 
AM
.,
2021
.
Model reduction and neural networks for parametric PDEs
.
SMAI J Comput Math
,
7
:
121
57
.

Carcione
 
JM
.,
1993
.
Seismic modeling in viscoelastic media
.
Geophysics
,
58
:
110
20
.

Carcione
 
JM
,
Herman
 
GC
,
Ten Kroode
 
APE
.,
2002
.
Seismic modelling
.
Geophysics
,
67
:
1304
25
.

Chen
 
H
,
Ge
 
Z
.,
2023
.
Physical information neural networks for 2D and 3D nonlinear Biot model and simulation on the pressure of brain
.
J Comput Phys
,
490
:
112309
.

Chen
 
H
,
Zhou
 
H
,
Jiang
 
S
,
Rao
 
Y
.,
2019
.
Fractional Laplacian viscoacoustic wave equation low-rank temporal extrapolation
.
IEEE Access
,
7
:
93187
97
.

Duveneck
 
E
,
Bakker
 
PM
.,
2011
.
Stable P-wave modeling for reverse-time migration in tilted TI media
.
Geophysics
,
76
:
S65
75
.

Grady
 
TJ
,
Khan
 
R
,
Louboutin
 
M
,
Yin
 
Z
,
Witte
 
PA
,
Chandra
 
R
,
Hewett
 
RJ
 et al. ,
2023
.
Model-parallel Fourier neural operators as learned surrogates for large-scale parametric PDEs
.
Comput Geosci
,
178
:
105402
.

Guo
 
P
,
McMechan
 
GA
,
Guan
 
H
.,
2016
.
Comparison of two viscoacoustic propagators for Q-compensated reverse time migration
.
Geophysics
,
81
:
S281
97
.

He
 
K
,
Zhang
 
X
,
Ren
 
S
,
Sun
 
J.
,
2016
.
Deep residual learning for image recognition
. In
Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)
, pp. 770–8.

Jiang
 
Z
,
Tahmasebi
 
P
,
Mao
 
Z
.,
2021
.
Deep residual U-net convolution neural networks with autoregressive strategy for fluid flow predictions in large-scale geosystems
.
Adv Water Resour
,
150
:
103878
.

Kalyani
 
VK
,
Pallavika
,
Chakraborty
 
SK
.,
2014
.
Finite-difference time-domain method for modelling of seismic wave propagation in viscoelastic media
.
Appl Math Comput
,
237
:
133
45
.

Karniadakis
 
GE
,
Kevrekidis
 
IG
,
Lu
 
L
,
Perdikaris
 
P
,
Wang
 
S
,
Yang
 
L
.,
2021
.
Physics-informed machine learning
.
Nat Rev Phys
,
3
:
422
40
.

Kosloff
 
DD
,
Baysal
 
E
.,
1982
.
Forward modeling by a Fourier method
.
Geophysics
,
47
:
1402
12
.

Krizhevsky
 
A
,
Sutskever
 
I
,
Hinton
 
GE
.,
2012
.
ImageNet classification with deep convolutional neural networks
.
Advances in Neural Information Processing Systems
,
25
.

LeCun
 
Y
,
Bengio
 
Y
.,
2015
.
Hinton
 
G
.
Deep learning
.
Nature
,
521
:
436
44
.

Li
 
Q
,
Zhou
 
H
,
Zhang
 
Q
,
Chen
 
H
,
Sheng
 
S
.,
2016
.
Efficient reverse time migration based on fractional Laplacian viscoacoustic wave equation
.
Geophys J Int
,
204
:
488
504
.

Li
 
Z
,
Kovachki
 
N
,
Azizzadenesheli
 
K
 et al. ,
2020b
.
Multipole graph neural operator for parametric partial differential equations
.
Advances in Neural Information Processing Systems
,
33
: 6755–66.

Li
 
Z
,
Kovachki
 
N
,
Azizzadenesheli
 
K
,
Liu
 
B
,
Bhattacharya
 
K
,
Stuart
 
A
,
Anandkumar
 
A.
,
2020a
.
Fourier neural operator for parametric partial differential equations
.
arXiv
.

Liu
 
Y
,
Sen
 
MK
.,
2009a
.
A new time–space domain high-order finite-difference method for the acoustic wave equation
.
J Comput Phys
,
228
:
8779
806
.

Liu
 
Y
,
Sen
 
MK
.,
2009b
.
A practical implicit finite-difference method: examples from seismic modelling
.
J Geophys Eng
,
6
:
231
49
.

Liu
 
Y
,
Sen
 
MK
.,
2011
.
Finite-difference modeling with adaptive variable-length spatial operators
.
Geophysics
,
76
:
T79
89
.

Long
 
J
,
Shelhamer
 
E
,
Darrell
 
T
.,
2015
.
Fully convolutional networks for semantic segmentation
. In
Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)
, pp. 3431–40.

Lu
 
L
,
Jin
 
P
,
Karniadakis
 
GE
.,
2019
.
Deeponet: learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators
.
arXiv
.

Moseley
 
B
,
Markham
 
A
,
Nissen-Meyer
 
T
.,
2020
.
Solving the wave equation with physics-informed deep learning
.
arXiv
.

Patnaik
 
S
,
Hollkamp
 
J P
.,
2020
.
Semperlotti
 
F
.
Applications of variable-order fractional operators: a review
.
Proc R Soc A
,
476
:
20190498
.

Qu
 
J
,
Cai
 
W
,
Zhao
 
Y
.,
2022
.
Learning time-dependent PDEs with a linear and nonlinear separate convolutional neural network
.
J Comput Phys
,
453
:
110928
.

Raissi
 
M
,
Perdikaris
 
P
,
Karniadakis
 
GE
.,
2019
.
Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations
.
J Comput Phys
,
378
:
686
707
.

Rashid
 
MM
,
Pittie
 
T
,
Chakraborty
 
S
 et al. ,
2022
.
Learning the stress-strain fields in digital composites using Fourier neural operator
.
Iscience
,
25
:
150452
.

Ren
 
Z
,
Liu
 
Y
.,
2013
.
A hybrid absorbing boundary condition for frequency-domain finite-difference modelling
.
J Geophys Eng
,
10
:
054003
.

Robertsson
 
JOA
,
Blanch
 
JO
,
Symes
 
WW
.,
1994
.
Viscoelastic finite-difference modelling
.
Geophysics
,
59
:
1444
56
.

Smith
 
JD
,
Azizzadenesheli
 
K
,
Ross
 
ZE
.,
2020
.
Eikonet: solving the eikonal equation with deep neural networks
.
IEEE Trans Geosci Remote Sens
,
59
:
10685
96
.

Song
 
C
,
Alkhalifah
 
T
,
Waheed
 
UB
.,
2021
.
Solving the frequency-domain acoustic VTI wave equation using physics-informed neural networks
.
Geophys J Int
,
225
:
846
59
.

Song
 
C
,
Wang
 
Y
.,
2022
.
High-frequency wavefield extrapolation using the Fourier neural operator
.
J Geophys Eng
,
19
:
269
82
.

Stoffa
 
PL
,
Fokkema
 
JT
,
de Luna Freire
 
RM
 et al. ,
1990
.
Split-step Fourier migration
.
Geophysics
,
55
:
410
21
.

Sun
 
J
,
Zhu
 
T
,
Fomel
 
S
.,
2015
.
Viscoacoustic modeling and imaging using low-rank approximation
.
Geophysics
,
80
:
A103
8
.

Wang
 
JX
,
Wu
 
JL
,
Xiao
 
H
.,
2017
.
Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data
.
Phys Rev Fluids
,
2
:
034603
.

Wang
 
N
,
Xing
 
G
,
Zhu
 
T
,
Zhou
 
H
,
Shi
 
Y
.,
2022
.
Propagating seismic waves in VTI attenuating media using fractional viscoelastic wave equation
.
J Geophys Res Solid Earth
,
127
:
e2021JB023280
.

Wang
 
N
,
Zhou
 
H
,
Chen
 
H
,
Xia
 
M
,
Wang
 
S
,
Fang
 
J
,
Sun
 
P
.,
2018
.
A constant fractional-order viscoelastic wave equation and its numerical simulation scheme
.
Geophysics
,
83
:
T39
48
.

Wei
 
W
,
Fu
 
LY
.,
2022
.
Small-data-driven fast seismic simulations for complex media using physics-informed Fourier neural operators
.
Geophysics
,
87
:
T435
46
.

Wen
 
G
,
Li
 
Z
,
Azizzadenesheli
 
K
,
Anandkumar
 
A
,
Benson
 
SM
.,
2022
.
U-FNO—an enhanced Fourier neural operator-based deep-learning model for multiphase flow
.
Adv Water Resour
,
163
:
104180
.

Wen
 
G
,
Tang
 
M
,
Benson
 
SM
.,
2021
.
Towards a predictor for CO2 plume migration using deep neural networks
.
Int J Greenh Gas Control
,
105
:
103223
.

Xing
 
G
,
Zhu
 
T
.,
2019
.
Modeling frequency-independent Q viscoacoustic wave propagation in heterogeneous media
.
J Geophys Res Solid Earth
,
124
:
11568
84
.

Xiong
 
Y
,
Guo
 
X
.,
2022
.
A short-memory operator splitting scheme for constant-Q viscoelastic wave equation
.
J Comput Phys
,
449
:
110796
.

Xu
 
J
,
Hu
 
H
,
Liu
 
QH
,
Han
 
B
.,
2023
.
A high-order perfectly matched layer scheme for second-order spectral-element time-domain elastic wave modelling
.
J Comput Phys
,
491
:
112373
.

Yang
 
Y
,
Gao
 
AF
,
Azizzadenesheli
 
K
,
Clayton
 
RW
,
Ross
 
ZE
.
2023
.
Rapid seismic waveform modeling and inversion with neural operators
.
IEEE Trans Geosci Remote Sens
,
61
:
1
12
.

Yao
 
J
,
Zhu
 
T
,
Hussain
 
F
,
Kouri
 
DJ
.,
2017
.
Locally solving fractional Laplacian viscoacoustic wave equation using Hermite distributed approximating functional method
.
Geophysics
,
82
:
T59
67
.

Yin
 
Z
,
Orozco
 
R
,
Louboutin
 
M
,
Herrmann
 
FJ
.,
2023
.
Solving multiphysics-based inverse problems with learned surrogates and constraints
.
Adv Model Simul Eng Sci
,
10
:
14
.

Zhang
 
T
,
Trad
 
D
,
Innanen
 
K
.,
2023
.
Learning to solve the elastic wave equation with Fourier neural operators
.
Geophysics
,
88
:
T101
19
.

Zhang
 
Y
,
Liu
 
Y
,
Xu
 
S
.,
2020
.
Arbitrary-order Taylor series expansion-based viscoacoustic wavefield simulation in 3D vertical transversely isotropic media
.
Geophys Prospect
,
68
:
2379
99
.

Zhang
 
Y
,
Liu
 
Y
,
Zhu
 
H
,
Chen
 
T
,
Li
 
J
.,
2023
.
Modified viscoelastic wavefield simulations in the time domain using the new fractional Laplacians
.
J Geophys Eng
,
19
:
346
61
.

Zhong
 
Z
,
Ni
 
B
,
Shi
 
Y
,
Shen
 
J
,
Du
 
X
.,
2023
.
Convolutional neural network-based seismic fragility analysis of subway station structure considering spatial variation of site shear-wave velocity
.
Comput Geotech
,
163
:
105741
.

Zhou
 
H
,
Liu
 
Y
,
Wang
 
J
.,
2023
.
K-space dispersion error compensators for the fractional spatial derivatives based constant-Q viscoelastic wave equation modelling
.
J Comput Phys
,
487
:
112161
.

Zhu
 
T
.,
2017
.
Numerical simulation of seismic wave propagation in viscoelastic-anisotropic media using frequency-independent Q wave equation
.
Geophysics
,
82
:
WA1
WA10
.

Zhu
 
T
,
Carcione
 
JM
.,
2014
.
Theory and modelling of constant-Q P- and S-waves using fractional spatial derivatives
.
Geophys J Int
,
196
:
1787
95
.

Zhu
 
Y
,
Zabaras
 
N
,
Koutsourelakis
 
P-S
,
Perdikaris
 
P
.,
2019
.
Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data
.
J Comput Phys
,
394
:
56
81
.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.