3D Acoustic Wave Equation Research Articles

A stable decoupled perfectly matched layer for the 3D wave equation using the nodal discontinuous Galerkin method

In outdoor acoustics, the calculations of sound propagating in air can be computationally heavy if the domain is chosen large enough for the sound waves to decay. The computational cost is lowered by strategically truncating the computational domain with an efficient boundary treatment. One commonly used boundary treatment is the perfectly matched layer (PML), which dampens outgoing waves without polluting the computed solution in the inner domain. The purpose of this study is to propose and assess a new perfectly matched layer formulation for the 3D acoustic wave equation, using the nodal discontinuous Galerkin finite element method. The formulation is based on an efficient PML formulation that can be decoupled to further increase the computational efficiency and guarantee stability without sacrificing accuracy. This decoupled PML formulation is demonstrated to be long-time stable, and an optimization procedure for the damping functions is proposed to enhance the performance of the formulation.

A robust numerical method based on a deep-learning operator for solving the 2D acoustic wave equation

The numerical solution of a wave equation plays a crucial role in computational geophysics problems, which forms the foundation of inverse problems and directly impacts the high-precision imaging results of earth models. However, common numerical methods often lead to significant computational and storage requirements. Due to the heavy reliance on forward-modeling methods in inversion techniques, particularly full-waveform inversion, enhancing the computational efficiency and reducing the storage demands of traditional numerical methods have become key issues in computational geophysics. We develop the deep Lax-Wendroff correction (DeepLWC) method, a deep-learning-based numerical format for solving 2D hyperbolic wave equations. DeepLWC combines the advantages of traditional numerical schemes with a deep neural network. We provide a detailed comparison of this method with the representative traditional Lax-Wendroff correction method. Our numerical results indicate that the DeepLWC significantly improves the calculation speed (by more than 10 times) and reduces the space needed for storage by more than 10,000 times compared with traditional numerical methods. In contrast to the more popular physics-informed neural network method, DeepLWC maximizes the advantages of traditional mathematical methods in solving partial differential equations and uses a new sampling approach, leading to improved accuracy and faster computations. It is particularly worth pointing out that DeepLWC introduces a novel research paradigm for numerical equation solving, which can be combined with various traditional numerical methods, enabling acceleration and reduction in the storage requirements of conventional approaches.

Robust data driven discovery of a seismic wave equation

SUMMARY Despite the fact that our physical observations can often be described by derived physical laws, such as the wave equation, in many cases, we observe data that do not match the laws or have not been described physically yet. Therefore recently, a branch of machine learning has been devoted to the discovery of physical laws from data. We test this approach for discovering the wave equation from the observed spatial-temporal wavefields. The algorithm first pre-trains a neural network (NN) in a supervised fashion to establish the mapping between the spatial-temporal locations (x, y, z, t) and the observation displacement wavefield function u(x, y, z, t). The trained NN serves to generate metadata and provide the time and spatial derivatives of the wavefield (e.g. utt and uxx) by automatic differentiation. Then, a preliminary library of potential terms for the wave equation is optimized from an overcomplete library by using a genetic algorithm. We, then, use a physics-informed information criterion to evaluate the precision and parsimony of potential equations in the preliminary library and determine the best structure of the wave equation. Finally, we train the ‘physics-informed’ neural network to identify the corresponding coefficients of each functional term. Examples in discovering the 2-D acoustic wave equation validate the feasibility and effectiveness of our implementation. We also verify the robustness of this method by testing it on noisy and sparsely acquired wavefield data.

A Splitting Nearly-Analytic Symplectic Partitioned Runge–Kutta Method for Solving 2D Elastic Wave Equations

In this paper, we present a new splitting nearly-analytic symplectic partitioned Runge–Kutta (SNSPRK) method for the two-dimensional (2D) elastic wave equations. It is an extension to elastic wave equation of our recent work on the locally one-dimensional nearly-analytic symplectic partitioned Runge–Kutta (LOD-NSPRK) method for the 2D acoustic wave equations. The method is based on the spatial differential operator-split technique, in which the resulting spatial discrete matrices are symmetric unlike the conventional nearly-analytic symplectic partitioned Runge–Kutta (NSPRK) method. The stability condition is given, which has more relax restriction for the time step than the NSPRK schemes. To show the performance of the new method, numerical experiments are given. Numerical results illustrate that the SNSPRK method has better long-term calculation capability as time proceeds than the NSPRK method.

Simulation of Wave Propagation Using Finite Differences in Oil Exploration

This paper presents a numerical solution for the 2D acoustic wave equation, considering heterogeneous media. It has been developed through a software in Fortran 90 that uses a second-order finite difference approximation. This program generates a set of patterns to detect the presence of oil in the subsurface. The algorithm is based on a geological domain where the sources (shots) and receivers are located. Each process takes care of a subset of sources and returns to the primary method patterns and seismograms corresponding to its group of sources. In the end, an image of the resulting seismogram is shown along the analyzed geologic profile. Stability and convergence tests were performed to ensure the reliability of the results. These tests were performed using a geological profile 100,000 m long and 17,400 m deep, divided into strata. For the execution of the software, a cluster of 16 processors was used as a computational platform.

Physics-informed neural networks for transcranial ultrasound wave propagation

Transcranial ultrasound imaging has been playing an increasingly important role in the non-invasive treatment of brain disorders. However, the conventional mesh-based numerical wave solvers, which are an integral part of imaging algorithms, suffer from limitations such as high computational cost and discretization error in predicting the wavefield passing through the skull. In this paper, we explore the use of physics-informed neural networks (PINNs) for predicting the transcranial ultrasound wave propagation. The wave equation, two sets of time snapshots data and a boundary condition (BC) are embedded as physical constraints in the loss function during training. The proposed approach has been validated by solving the two-dimensional (2D) acoustic wave equation under three increasingly complex spatially varying velocity models. Our cases demonstrate that due to the meshless nature of PINNs, they can be flexibly applied to different wave equations and types of BCs. By adding physics constraints to the loss function, PINNs can predict wavefields far outside the training data, providing ideas for improving the generalization capability of existing deep learning methods. The proposed approach offers exciting perspectives because of the powerful framework and simple implementation. We conclude with a summary of the strengths, limitations and further research directions of this work.

A Compact High-Order Finite-Difference Method with Optimized Coefficients for 2D Acoustic Wave Equation

High-precision finite difference (FD) wavefield simulation is one of the key steps for the successful implementation of full-waveform inversion and reverse time migration. Most explicit FD schemes for solving seismic wave equations are not compact, which leads to difficulty and low efficiency in boundary condition treatment. Firstly, we review a family of tridiagonal compact FD (CFD) schemes of various orders and derive the corresponding optimization schemes by minimizing the error between the true and numerical wavenumber. Then, the optimized CFD (OCFD) schemes and a second-order central FD scheme are used to approximate the spatial and temporal derivatives of the 2D acoustic wave equation, respectively. The accuracy curves display that the CFD schemes are superior to the central FD schemes of the same order, and the OCFD schemes outperform the CFD schemes in certain wavenumber ranges. The dispersion analysis and a homogeneous model test indicate that increasing the upper limit of the integral function helps to reduce the spatial error but is not conducive to ensuring temporal accuracy. Furthermore, we examine the accuracy of the OCFD schemes in the wavefield modeling of complex structures using a Marmousi model. The results demonstrate that the OCFD4 schemes are capable of providing a more accurate wavefield than the CFD4 scheme when the upper limit of the integral function is 0.5π and 0.75π.

Modeling seismic waves in ocean with the presence of irregular seabed and rough sea surface

Abstract For the simulation of seismic wavefields in the ocean, it is very important to consider the influences of the irregular seabed and the rough sea surface. To numerically model seismic waves in the ocean, the fluid seawater and the solid seabed must be treated separately, and the air–fluid and fluid–solid boundary conditions must be considered simultaneously. A fluid–solid configuration is created with the influences of a rough sea surface and an irregular seabed to carry out wavefield simulation. Based on the boundary conformal grids, the 2D acoustic wave equations (in the fluid media), the 2D elastic wave equations (in the solid media) and two boundary conditions (at the interface between air and fluid and the interface between fluid and solid) are converted from the Cartesian coordinate system into a curvilinear coordinate system. In curvilinear coordinates, the seismic wave equations and boundary conditions are discretized with finite-difference operators. At the left, right and bottom boundaries of the model, the convolutional perfectly matched layer boundary condition is used to absorb artificial reflections. The numerical simulation results show that the snapshots and shot gathers of the rough sea surface model are more complex and contain stronger reflections than those of the flat model. The rough sea surface and the pattern of the irregular seabed have a great influence on the propagation of seismic waves in the ocean.

An optimized three-dimensional time-space domain staggered-grid finite-difference method

Numerical simulation of three-dimensional (3D) seismic wavefields forms the basis of the research on the migration methods of 3D seismic data based on wave equations. Because the simulation precision of wavefield extrapolation determines the imaging accuracy to a certain extent, it is very important to study how to enhance the forward modeling precision of 3D seismic wavefields. Thus, we build on an optimized 3D staggered-grid finite-difference (SFD) method with high simulation precision based on two-dimensional (2D) seismic modeling. Since it generates the corresponding difference coefficients by utilizing the least square (LS) method to minimize the objective function constructed by the time-space domain dispersion relation of the 3D acoustic wave equation, our optimized time-space domain LS-based 3D SFD method can effectively enhance the modeling precision of the 3D seismic wavefields in theory compared with the 3D SFD methods based on the Taylor-series expansion (TE), especially for the large wavenumber range. Examining the numerical dispersion, algorithm stability and computational cost, we compare our optimized time-space domain LS-based 3D SFD method with three conventional TE-based and LS-based 3D SFD methods to illustrate and demonstrate its effectiveness and feasibility. The numerical examples from different 3D models suggest that our optimized time-space domain LS-based 3D SFD method can generate less numerical dispersion and higher simulation accuracy for 3D seismic wavefields than three other conventional 3D SFD methods, but its stability condition is stricter and its computational cost is slightly higher.

An accurate and efficient local one-dimensional method for the 3D acoustic wave equation

Abstract We establish an accurate and efficient scheme with four-order accuracy for solving three-dimensional (3D) acoustic wave equation. First, the local one-dimensional method is used to transfer the 3D wave equation into three one-dimensional wave equations. Then, a new scheme is obtained by the Padé formulas for computation of spatial second derivatives and the correction of the truncation error remainder for discretization of temporal second derivative. It is compact and can be solved directly by the Thomas algorithm. Subsequently, the Fourier analysis method and the Lax equivalence theorem are employed to prove the stability and convergence of the present scheme, which shows that it is conditionally stable and convergent, and the stability condition is superior to that of most existing numerical methods of equivalent order of accuracy in the literature. It allows us to reduce computational cost with relatively large time step lengths. Finally, numerical examples have demonstrated high accuracy, stability, and efficiency of our method.

A High Accuracy Local One-Dimensional Explicit Compact Scheme for the 2D Acoustic Wave Equation

In this paper, we develop a highly accurate and efficient finite difference scheme for solving the two-dimensional (2D) wave equation. Based on the local one-dimensional (LOD) method and Padé difference approximation, a fourth-order accuracy explicit compact difference scheme is proposed. Then, the Fourier analysis method is used to analyze the stability of the scheme, which shows that the new scheme is conditionally stable and the Courant-Friedrichs-Lewy (CFL) condition is superior to most existing methods of equivalent order of accuracy in the literature. Finally, numerical experiments demonstrate the high accuracy, stability, and efficiency of the proposed method.

Upscaling acoustic wave equation using renormalization group theory

Seismic waves interact with a broad range of heterogeneities as they propagate through the earth. Simulating this full range of scales for wave propagation requires capturing heterogeneities of all scales, which can be computationally unaffordable. In such cases, we have relied on macroscopic representations of media obtained through an upscaling process that preserves the effects of small-scale heterogeneities (in comparison with the wavelengths of interest). Here, we discuss the application of the renormalization group (RG) theory-based upscaling to the 2D acoustic wave equation. RG-based upscaling requires constructing a special Fourier operator and is implemented using a domain decomposition in conjunction with an “expansion and truncation” method to mitigate edge effects. We test this upscaling method on several benchmark models and in the context of reverse time migration. The upscaled models obtained using this method indicate a good consistency for generated waveforms, whereas the runtime for simulations is reduced by at least an order of magnitude.

2D acoustic wavefield simulation by improved stiffness and mass matrices of bilinear square element

Abstract In the frequency domain, we optimize the entries of the element stiffness matrix computed on a square element with the purpose of approximating the 2D acoustic wave equation with the second spatial order. The optimized matrices are computed from the minimization of the normalized phase velocity as a function of propagation angle and the number of grid points per wavelength, leading to a remarkable reduction in the numerical error. The optimized number of 4.1 sample points in a wavelength, with a maximum error in velocities <0.3%, offering the simulation of large-scale realistic models. The superior performance of the optimized matrices to those from the lumped, consistent and eclectic matrices, is evident from the numerical examples presented, with consequent reduction not only in the numerical dispersion/anisotropy but also in an improved resolution. It is noteworthy that optimized matrices give modeling results with better accuracy than models from the high-order spectral element method. Our methodology is easily extendable to the accurate discretizing of the general 3D elastic wave equation that optimizes the bandwidth of the impedance matrix, availing computational resources for accurately modeling the compressional and shear wave velocities, density, as well as multi-parameter inversion.

Optimized Mass and Stiffness Matrices for Accurate and Efficient 3D Acoustic Wave Modeling

AbstractNumerical modeling of the 3D acoustic wave equation regularly faces a trade‐off between accuracy and efficiency. Usually, an improved computational efficiency leads to a compromised numerical accuracy, in which the numerical dispersion error increases. Such a trade‐off is commonly reported for the finite difference methods, and the consistent, lumped, eclectic finite element matrices, and the high‐order spectral element methods. Therefore, we present a methodology that simultaneously optimizes both the accuracy and the computational efficiency of the finite element method for modeling the 3D acoustic wave equation in the frequency domain. The numerical dispersion objective function is derived and minimized for a tri‐linear interpolation function of a cube element in terms of the mass and the stiffness matrices. With four grid points per wavelength sampling, for less than 0.5% error in velocities, the optimized modeled solutions outperform the results from extant approaches. Numerical experiments with the homogeneous, 3D unbounded 2‐layer and the 3D Society of Exploration Geophysicists salt velocity models reveal that the optimized accuracy is persistent across all source‐receiver offsets, a performance observably absent from other finite element matrices. The optimized matrices are not only useful for numerical modeling, but are also readily available for accurate computation of the sensitivity matrices of the density and the bulk modulus parameters for the multi‐parameterization inversion.

Seismic Wave Propagation and Inversion with Neural Operators

AbstractSeismic wave propagation forms the basis for most aspects of seismological research, yet solving the wave equation is a major computational burden that inhibits the progress of research. This is exacerbated by the fact that new simulations must be performed whenever the velocity structure or source location is perturbed. Here, we explore a prototype framework for learning general solutions using a recently developed machine learning paradigm called neural operator. A trained neural operator can compute a solution in negligible time for any velocity structure or source location. We develop a scheme to train neural operators on an ensemble of simulations performed with random velocity models and source locations. As neural operators are grid free, it is possible to evaluate solutions on higher resolution velocity models than trained on, providing additional computational efficiency. We illustrate the method with the 2D acoustic wave equation and demonstrate the method’s applicability to seismic tomography, using reverse-mode automatic differentiation to compute gradients of the wavefield with respect to the velocity structure. The developed procedure is nearly an order of magnitude faster than using conventional numerical methods for full waveform inversion.

Adaptive 27-point finite-difference frequency-domain method for wave simulation of 3D acoustic wave equation

The finite-difference frequency-domain (FDFD) method has important applications in the wave simulation of various wave equations. To promote the accuracy and efficiency for wave simulation with the FDFD method, we have developed a new 27-point FDFD stencil for 3D acoustic wave equations. In the developed stencil, the FDFD coefficients depend not only on the ratios of cell sizes in the x-, y-, and z-directions but also depend on the spatial sampling density (SD) in terms of the number of wavelengths per grid. The corresponding FDFD coefficients can be determined efficiently by making use of the plane-wave expression and the lookup table technique. We also develop a new way for designing an adaptive FDFD stencil by directly adding some correction terms to the conventional second-order FDFD stencil, which is simpler to use and easier to generalize. Corresponding dispersion analysis indicates that, compared to the optimal 27-point stencil derived from the average-derivative method (ADM), the developed adaptive 27-point stencil can reduce the required SD from approximately 4 to 2.2 points per wavelength (PPW) for a cubic mesh and to 2.7 PPW for a general cuboid mesh. Numerical examples of a 3D homogeneous model and SEG/EAGE salt-dome model indicate that the developed stencil is more accurate than the ADM 27-point stencil for cubic and general cuboid meshes, while requiring similar CPU time and computational memory as the ADM 27-point stencil for direct and iterative solvers.

Time-space domain dispersion reduction schemes in the uniform norm for the 2D acoustic wave equation

Finite-difference (FD) methods for the wave equation are flexible, robust and easy to implement. However, they in general suffer from numerical dispersion. FD methods based on accuracy give good dispersion at low frequencies, but waves tend to disperse for higher wavenumbers. Moreover, waves in higher dimensions also suffer from dispersion errors in all propagation angles. In this work, we give a unified methodology to derive dispersion reduction FD schemes for the two dimensional acoustic wave equation. This new methodology would generate schemes that give the theoretical minimum dispersion error in the uniform norm. Stability criteria are discussed for the general scheme. We also motivate the equivalence of popular dispersion reduction methodology and theoretical numerical reduction methodology.

Applying an advanced temporal and spatial high-order finite-difference stencil to 3D seismic wave modeling

Present finite-difference (FD) algorithms for modeling seismic wave propagation in 3D acoustic media are mainly based on the traditional orthogonal stencil and can achieve spatial high-order but temporal second-order accuracy. Therefore, these approaches may result in visible temporal dispersion and even instability for relatively large time sampling intervals. In this paper, we develop an advanced FD stencil with temporal and spatial arbitrary even-order accuracy to solve the 3D acoustic wave equation. The temporal derivative is approximated by an octahedron stencil with the operator length N together with the traditional second-order FD. Meanwhile, the spatial derivatives are calculated by the orthogonal stencil with the operator length M. The plane-wave theory is introduced to derive the time-space-domain dispersion relation and estimate the corresponding high-order FD coefficients by Taylor expansion. To achieve higher accuracy, we further optimize FD coefficients by applying least-squares (LS) to minimize the relative error of the time-space-domain dispersion relation. Dispersion and stability analyses show that our new schemes have higher modeling accuracy and propagation stability than the conventional method. Numerical examples demonstrate that the proposed FDs can generate greater temporal accuracy on the prerequisite of preserving satisfactory spatial accuracy. Our new temporal high-order FD stencil by incorporating larger time step, shorter operator lengths, LS-based coefficients and variable-length operators can be regarded as a kind of accurate and efficient tool for seismic wave modeling.

Optimizing orthogonal-octahedron finite-difference scheme for 3D acoustic wave modeling by combination of Taylor-series expansion and Remez exchange method

Using orthogonal-octahedron finite-difference (FD) stencils, even-order accuracy for temporal and spatial derivatives can be reached for 3D acoustic wave equation modelling. However, the required length of FD operator is still long for high accuracy modelling. To tackle this issue, we have developed an optimization method using a combination of Taylor-series expansion plus Remez exchange methods. The Taylor-series expansion for the octahedron-shaped operator ensures high temporal accuracy and the Remez exchange for axial operator further improves the spatial accuracy. The comparison between our method, exisiting 3D temporal high order method and space domain least-square opimization method validates the efficiency and accuracy superiority of the new method.

A parallel hybrid implementation of the 2D acoustic wave equation

Abstract In this paper, we propose a hybrid parallel programming approach for a numerical solution of a two-dimensional acoustic wave equation using an implicit difference scheme for a single computer. The calculations are carried out in an implicit finite difference scheme. First, we transform the differential equation into an implicit finite-difference equation and then using the alternating direction implicit (ADI) method, we split the equation into two sub-equations. Using the cyclic reduction algorithm, we calculate an approximate solution. Finally, we change this algorithm to parallelize on graphics processing unit (GPU), GPU + Open Multi-Processing (OpenMP), and Hybrid (GPU + OpenMP + message passing interface (MPI)) computing platforms. The special focus is on improving the performance of the parallel algorithms to calculate the acceleration based on the execution time. We show that the code that runs on the hybrid approach gives the expected results by comparing our results to those obtained by running the same simulation on a classical processor core, Compute Unified Device Architecture (CUDA), and CUDA + OpenMP implementations.