Effective filtering of artifacts for finite‐difference implementations of paraxial wave equation migration

Heterogeneous wave equations are more complicated numerically than homogeneous wave equations, but are necessary for physical validity. A wide variety of numerical solutions of seismic wave equations is available, but most produce strong numerical artefacts and local instabilities where model parameters change rapidly. Accuracy and stability of heterogeneous equations is achieved through staggered‐grid formulations. A new pseudospectral staggered‐grid algorithm is developed for the poroelastic (Biot) equations. The algorithm may be reduced to handle the elastic and acoustic limits of the Biot equations. Comparisons of results from poroelastic, elastic, acoustic and scalar computations for a 2D model show that porous medium parameters may affect amplitudes significantly. The use of homogeneous wave equations for modelling of a heterogeneous medium, or of a centred rather than a staggered grid, or of simplified (e.g. acoustic) wave equations when elastic or poroelastic media are synthesized, may produce erroneous or ambiguous interpretations.

We study monochromatic, scalar solutions of the Helmholtz and paraxial wave equations (PWEs) from a field-theoretic point of view. We introduce appropriate time-independent Lagrangian densities for which the Euler–Lagrange equations reproduces either Helmholtz and PWEs with the z-coordinate, associated with the main direction of propagation of the fields, playing the same role of time in standard Lagrangian theory. For both Helmholtz and paraxial scalar fields, we calculate the canonical energy-momentum tensor and determine the continuity equations relating ‘energy’ and ‘momentum’ of the fields. Eventually, the reduction of the Helmholtz wave equation to a useful first-order Dirac form, is presented. This work sheds some light on the intriguing and not so acknowledged connections between angular spectrum representation of optical wavefields, cosmological models and physics of black holes.

The numerlcal sUmulation of the elastic wave equation has been become one of the fundamental seismic modeling methods. In recent years, because high speed and large storage computers were avaiable, the numerical stimulation of the Blot’s coupled equation was developed. Biot’s equation takes account of energy dissipation due to the relative motion between a viscous pore fluid and the porous solid matrix. But for the samc size model the numerical s!hnulaUon of Biot’s equation takes much more computer storage and hours than that of elastic wave equation. The realistic oilfields are composed of non-porous rocks (elastic media) and fluid-saturated porous rocks (Blot’s media). In this paper, we pmaent a hybrid finite difference scheme to stimulate wave propagation in reallstic oilfields, which describes the wave propagation in an elastic medium by elastic wave equation, and in a fluid-saturated porous medium by Biot’s coupled wave equation. Numerical example is presented to demonstrate the scheme. By comparing the result with that computed by acoustic wave equation, the effects of parameters uch as porosity and fluid content can be investigated.

SUMMARY We first transform the seismic wave equations in 2-D inhomogeneous anisotropic media into a system of first-order partial differential equations with respect to time t. And then we develop a new explicit split-step algorithm (SSA) to solve the transformed equations through using the implicit Adams method and high-order interpolation approximation. Our algorithm enables wave propagation to be simulated in two dimensions through generally isotropic and anisotropic models. The high-order space derivatives in this SSA are determined by using the wave displacement, particle-velocity and its gradients simultaneously, while the time advancing of modelling seismic propagation uses the slit-operator algorithm based on the third-order implicit Adams method. On the basis of such a structure, the so-called SSA can suppress effectively the numerical dispersion and source-generated noises caused by discretizing the wave equations when too-coarse grids are used, and is third-order accuracy in time and fourth-order accuracy in space. We explore the theoretical properties of the SSA including the stability criteria of the SSA for solving 1-D and 2-D scalar wave equations, numerical dispersion, theoretical error and numerical error and computational efficiency when using the SSA to model the acoustic wave fields. Meanwhile, we also present the synthetic seismograms generated by the SSA in the three-layer isotropic medium and the 2-D elastic two-layer isotropic and anisotropic wavefield snapshots computed by the fourth-order Lax–Wendroff correction (LWC) scheme and the SSA, respectively. Numerical calculations of the relative errors show that the numerical error of the SSA is less than those of the conventional finite difference (FD) method and the fourth-order LWC scheme. Promising numerical results further illustrate that the SSA has less numerical dispersions and can suppress effectively the source-generated noises, further resulting in both increasing the computational efficiency and saving the space storage.

We present a quantum algorithmic framework for simulating linear, anti-Hermitian (lossless) wave equations in heterogeneous, anisotropic, and time-independent media. This framework encompasses a broad class of wave equations, including the acoustic wave equation, Maxwell’s equations, and the elastic wave equation. Our formulation is compatible with standard numerical discretization schemes and allows for the efficient implementation of multiple practically relevant time- and space-dependent sources. Furthermore, we demonstrate that subspace energies can be extracted and wave fields compared through an l2 loss function, achieving optimal precision scaling with the number of samples taken. Additionally, we introduce techniques for incorporating boundary conditions and linear constraints that preserve the anti-Hermitian nature of the equations. Leveraging the Hamiltonian simulation algorithm, our framework achieves a quartic speedup over classical solvers in three-dimensional simulations, under conditions of sufficiently global measurements and compactly supported sources and initial conditions. This quartic speedup is optimal for time-domain solutions, as the Hamiltonian of the discretized wave equations has local couplings. In summary, our framework provides a versatile approach for simulating wave equations on quantum computers, offering substantial speedups over state-of-the-art classical methods.

—Endeavors to simulate the 3-D seismic wave propagation effects in the 2-D environment have stimulated the development of 2.5-D simulation techniques. Liner (1991) has developed a damped 2.5-D acoustic wave equation using the 3-D Green's function for a constant density medium. This paper presents the finite-difference simulation and study of various numerical artifacts, using Liner's acoustic wave equation. Comparisons of second-order as well as fourth-order accurate 2.5-D and 2-D acoustic wave simulation results have been given. The snapshots at different times and zero off-set response of a salt dome model have been computed and strong reflections and diffraction phases have been reported. Also, the analysis of amplitude behaviors and computational time for Liner's (1991) and Williamson and Pratt's (1995) 2.5-D wave equations, with and without the last term (1/t 2),reveals that this last term is superfluous as suggested by Stockwell (1995). The effectiveness of Sponge transmissive as well as Clayton and Engquist (1980) absorbing boundary conditions has been studied. A modified Ricker wavelet, the second derivative of the convolution of Gaussian function and a polynomial window, has been used as a source. It has been found that the stability condition and the requirement of the number of grid points per wavelength to avoid grid dispersion are the same as for the 2-D acoustic cases.

The transport-of-intensity equation (TIE) is often used to determine the phase and amplitude profile of a complex object by monitoring the intensities at different distances of propagation or around the image plane. TIE results from the imaginary part of the paraxial wave equation and is equivalent to the conservation of energy. The real part of the paraxial wave equation gives the eikonal equation in the presence of diffraction. Since propagation of the optical field between different planes is governed by the (paraxial) wave equation, both real and imaginary parts need to be satisfied at every propagation plane. In this work, the solution of the TIE is optimized by using the real part of the paraxial wave equation as a constraint. This technique is applied to the more exact determination of imaging the induced phase of a liquid heated by a focused laser beam, which has been previously computed using TIE only. Retrieval of imaged phase using the TIE is performed by using the constraint that naturally arises from the real part of the paraxial wave equation.

We construct exact solutions to the paraxial wave equation in free space characterized by stable caustics. First, we show that any solution of the paraxial wave equation can be written as the superposition of plane waves determined by both the Hamilton-Jacobi and Laplace equations in free space. Then using the five elementary stable catastrophes, we construct solutions of the Hamilton-Jacobi and Laplace equations, and the corresponding exact solutions of the paraxial wave equation. Therefore, the evolution of the intensity patterns is governed by the paraxial wave equation and that of the corresponding caustic by the Hamilton-Jacobi equation.

Full waveform inversion (FWI) is an optimization procedure which minimizes the waveform difference between the observed and calculated seismic data. The method, through an iterative forward modeling and inversion based on the acoustic or elastic wave equation, has recently gained wide interest with specific application in exploration and producing. FWI is wellknown for its ability to obtain high resolution and accurate subsurface models in areas with complex geological structures related to the presence of salt bodies encountered in deep water exploration. However, conventional FWI algorithms require seismic data that contains low frequencies (<5Hz). Lack of low frequencies results in unreliable velocity models for seismic imaging. In this paper, we aim to relax the low-frequency requirement of FWI and propose a technique which employs principles from sparse deconvolution in order to recover the missing low frequencies. Through a series of synthetic examples, it is shown that our proposed methodology is able to preserve and utilize the traveltime information of the original data; hence, it allows convergence of the FWI. The major advantage of this method is that it is fully automatic, involving no residual moveout or traveltime difference analysis.

In marine seismic exploration, ocean-bottom cable techniques accurately record the multicomponent seismic wavefield; however, the seismic wave propagation in fluid–solid media cannot be simulated by a single wave equation. In addition, when the seabed interface is irregular, traditional finite-difference schemes cannot simulate the seismic wave propagation across the irregular seabed interface. Therefore, an acoustic–elastic forward modeling and vector-based P- and S-wave separation method is proposed. In this method, we divide the fluid–solid elastic media with irregular interface into orthogonal grids and map the irregular interface in the Cartesian coordinates system into a horizontal interface in the curvilinear coordinates system of the computational domain using coordinates transformation. The acoustic and elastic wave equations in the curvilinear coordinates system are applied to the fluid and solid medium, respectively. At the irregular interface, the two equations are combined into an acoustic–elastic equation in the curvilinear coordinates system. We next introduce a full staggered-grid scheme to improve the stability of the numerical simulation. Thus, separate P- and S-wave equations in the curvilinear coordinates system are derived to realize the P- and S-wave separation method.

We have extended a new highly accurate numerical scheme for unstructured 2D and 3D meshes based on the discontinuous Galerkin approach to simulate seismic wave propagation in heterogeneous media containing fluid-solid interfaces. Because of the formulation of the wave equations as a unified first-order hyperbolic system in velocity stress, the fluid can be in movement along the interface. The governing equations within the moving fluid are derived from well-known first principles in fluid mechanics. The discontinuous Galerkin approach allows for jumps of the material parameters and the solution across element interfaces, which are handled by Riemann solvers or numerical fluxes. The use of Riemann solvers at the element interfaces makesthe treatment of the fluid particularly simple bysetting the shearmodulus in the fluid region to zero. No additional compatibility relations, such as vanishing shear stress or continuity of normal stresses, are necessary to couple the solid and fluid along an interface. The Riemann solver automatically recognizes the jump of the material coefficients at the interface and provides the correct numerical fluxes for fluid-solid contacts. Therefore, wave propagation in the entire computational domain containing heterogeneous media, namely moving fluids and elastic solids, can be described by a uniform set of acoustic and elastic wave equations. The accuracy of the proposed scheme is confirmed by comparing numerical results against analytic solutions. The potential of the new method was demonstrated in a 3D model problem typical for marine seismic exploration with a fluid-solid interface determined by a complicated bathymetry.

This paper addresses the direct and inverse source problems for the stochastic acoustic, biharmonic, electromagnetic, and elastic wave equations in a unified framework. The driven source is assumed to be a centered generalized microlocally isotropic Gaussian random field, whose covariance and relation operators are classical pseudodifferential operators. Given the random source, the direct problems are shown to be well-posed in the sense of distributions and the regularity of the solutions are given. For the inverse problems, we demonstrate by ergodicity that the principal symbols of the covariance and relation operators can be uniquely determined by a single realization of the far-field pattern averaged over the frequency band with probability one.

Forward seismic modeling is a useful tool for simulating the response of the low velocity zones near the ground surface and delicate structures at a deeper depth in northwestern China. The main objective of forward seismic modeling is to solve the wave equation accurately. The finite difference method (FDM) is a commonly used numerical method for solving the equation. It can obtain very high accuracy with high‐order difference scheme if the grid configuration is appropriate. However it shows numerical dispersion and its stability criteria is also very strict. In order to deal with these problems, a new precise integration method is advanced in this paper, using an integration scheme in the temporal domain rather than a traditional low‐order difference scheme. The precise integration method employs a difference approximation in the spatial domain, converting the acoustic or elastic wave equation into a set of differential equations with respect to time and then solves this equation system using an integration method. The final solution is attributed to an integration of a matrix of exponential functions. Theoretically this solution is accurate in the temporal domain, hence increasing the total accuracy. The whole spatial domain can be divided into several sub‐domains, considering computer capacity and speed and then the solutions put together. This revised sub‐domain precise integration method (SPIM) is more stable than the FDM and more accurate as is tested in the text. Besides it can greatly reduce the numerical dispersion with the same grid configuration as the FDM, though it may take a little more CPU time. A numerical test is presented to justify this method and then an application is shown to indicate that the SPIM is of potential value in the geophysical exploration field.

In recent years, half-precision floating-point arithmetic has gained wide support in hardware and software stack thanks to the advance of artificial intelligence and machine-learning applications. Operating at half precision can significantly reduce the memory footprint compared to operating at single or double precision. For memory-bound applications such as time-domain wave simulations, this is an attractive feature. However, the narrower width of the half-precision data format can lead to degradation of the solution quality due to larger roundoff errors. We illustrate the negative impact caused by the accumulation of roundoff errors in wave simulations with carefully designed numerical experiments. Specifically, the energy-conserving property of the wave equations is used as a convenient diagnosis tool. The corresponding remedy in the form of compensated sum is then provided, with its efficacy demonstrated using numerical examples with acoustic and elastic wave equations on hardware that support half-precision arithmetic natively.

Full-waveform inversion (FWI) is a non-destructive health monitoring technique that can be used to identify and quantify the embedded anomalies. The forward modeling of the FWI consists of a simulation of elastic wave equation to generate synthetic data. Thus the accuracy of the FWI method highly depends on the simulation method used in the forward modeling. Simulation of a 3-D seismic survey with small-scale heterogeneities is impossible with the classic finite difference approach even on modern super computers. In this work, we adopted a mesh refinement approach for simulation of the wave equation in the presence of small-scale heterogeneities. This approach uses cubic smoothing spline interpolation for spatial mesh refinement step in solving the wave equation. The simulation results for the 2-D elastic wave equation are presented and compared with the classic finite difference approach.