Search

The pseudo-differential operator with symbol | k |α has been widely used in seismic modelling and imaging when involving attenuation, anisotropy and one-way wave equation, which is usually calculated using the pseudo-spectral method. For large-scale problems, applying high-dimensional Fourier transforms to solve the wave equation that includes pseudo-differential operators is much more expensive than finite-difference approaches, and it is not suitable for parallel computing with domain decomposition. To mitigate this difficulty, we present a truncated space-domain convolution method to efficiently compute the pseudo-differential operator |$\sqrt{-\nabla ^2}$|⁠ , and then apply it to viscoacoustic reverse-time migration. Although |$\sqrt{-\nabla ^2}$| is theoretically non-local in the space domain, we take the limited frequency band of seismic data into account, and constrain the approximated convolution stencil to a finite length. The convolution coefficients are computed by solving a least-squares inverse problem in the wavenumber domain. In addition, we exploit the symmetry of the resulting convolution stencil and develop a fast spatial convolution algorithm. The applications of the proposed method in Q -compensated reverse-time migration demonstrate that it is a good alternative to the pseudo-spectral method for computing the pseudo-differential operator |$\sqrt{-\nabla ^2}$|⁠ , with almost the same accuracy but much higher efficiency. [ABSTRACT FROM AUTHOR]

By fitting observed data with predicted seismograms, least‐squares migration (LSM) computes a generalized inverse for a subsurface reflectivity model, which can improve image resolution and reduce artifacts caused by incomplete acquisition. However, the large computational cost of LSM required for simulations and migrations limits its wide applications for large‐scale imaging problems. Using point‐spread function (PSF) deconvolution, we present an efficient and stable high‐resolution imaging method. The PSFs are first computed on a coarse grid using local ray‐based Gaussian beam Born modeling and migration. Then, we interpolate the PSFs onto a fine‐image grid and apply a high‐dimensional Gaussian function to attenuate artifacts far away from the PSF centers. With 2D/3D partition of unity, we decompose the traditional adjoint migration results into local images with the same window size as the PSFs. Then, these local images are deconvolved by the PSFs in the wavenumber domain to reduce the effects of the band‐limited source function and compensate for irregular subsurface illumination. The final assembled image is obtained by applying the inverse of the partitions for the deconvolved local images. Numerical examples for both synthetic and field data demonstrate that the proposed PSF deconvolution can significantly improve image resolution and amplitudes for deep structures, while not being sensitive to velocity errors as the data‐domain LSM. Plain Language Summary: Seismic imaging is an important tool to detect hydrocarbon resource and study deep Earth's structure. Traditional ray‐based and wave equation imaging methods commonly extrapolate observed data and apply appropriate imaging condition to construct subsurface impedance interfaces. Mathematically, these methods can be considered as adjoint operators of seismic modeling, but they have difficulty to produce high‐quality images of complicated subsurface structures. To mitigate this issue, we present an efficient and stable high‐resolution imaging approach, i.e., the point‐spread function (PSF) deconvolution. The locally temporal and spatial limitations are utilized to efficiently compute the PSFs on a coarse grid using a ray‐based Gaussian beam propagator. Then, an on‐the‐fly interpolation on a fine grid and a wavenumber‐domain deconvolution are applied to improve image resolution. Numerical experiments for benchmark model and field data demonstrate the feasibility and adaptability of the proposed method for imaging complicated salt structure and low‐signal‐to‐noise data. Key Points: Conventional least‐squares migration is too expensive to be widely applied for large‐scale problems under current computational capacityAn efficient approach for calculating point‐spread functions using local modeling and migration is proposedThe point‐spread function deconvolution is then applied to incorporate the Hessian effect and improve image quality [ABSTRACT FROM AUTHOR]

Ray theory gives a high‐frequency asymptotic solution of the wave equation, which has been widely used in the seismic study because of its high computational efficiency and good physical insights for elementary waves. However, the fundamental high‐frequency assumption makes it difficult to accurately simulate the finite‐frequency effect of wave propagations. On the other hand, the Gaussian beam, as one of the advanced ray‐based methods, overcomes the amplitude singularity problem near the caustics by introducing a complex‐valued paraxial approximation. But Gaussian beams only use the velocity and up to second‐order velocity derivative of central rays and thus does not accurately simulate the response of heterogeneity away from the rays. To bridge the gap between the ray theory and wave‐equation methods, we present a two‐way beam wave method and apply it to seismic imaging. A fan of central rays are first shot to extract local ray‐centered model parameters in the global Cartesian coordinates. Then, a finite‐difference method is used to solve the acoustic wave equation in the ray‐centered coordinates to simulate wave propagations near the central rays. The beam wave method can accurately simulate the response of the velocity heterogeneities in the ray tubes and therefore produces more accurate wavefields. In addition, we discretize the ray‐centered wave equation onto a non‐orthogonal grid, which considerably reduces the computational cost. We apply the two‐way beam wave method in seismic imaging and develop a beam wave reverse‐time migration. It inherits the advantages of both ray‐based and wave equation migrations, and produces clear images for complicated structures with fewer artifacts than traditional imaging methods. Plain Language Summary: Seismic ray theory has evolved from classical asymptotic ray theory, through paraxial ray theory and Maslov's method, and then to the Gaussian beam method. Although the computational accuracy of the ray‐based methods for seismic modeling and imaging have been greatly improved, the fundamental high‐frequency approximation still prevents them from accurately describing finite‐frequency wave propagations in complex structures. We introduce the idea of reverse‐time migration into the beam method and develop a two‐way beam wave scheme for seismic modeling and imaging. Local velocities in the vicinity of central rays are first extracted and then the finite‐difference approach is used to solve the wave equation in the ray‐centered coordinate. Accurate wavefield simulation and beam‐to‐beam imaging condition enables us to achieve clear migrated images with fewer artifacts than traditional Gaussian beam migration and reverse‐time migration. Numerical examples for benchmark models and field data demonstrate the feasibility and adaptability of the proposed two‐way beam wave method. Key Points: It is difficult for ray theory to accurately simulate seismic wave phenomena in areas with strong velocity variationsA two‐way beam wave method is proposed to bridge the gap between ray‐based and wave equation approachesA beam wave reverse‐time migration approach is developed to improve the image quality for complicated subsurface environments [ABSTRACT FROM AUTHOR]

Accurately characterizing seismic attenuation effects on wave propagations is crucially important for structure interpretation and reservoir evaluation. The conventional fractional viscoelastic wave equation is not satisfactory on accuracy for small Q values. To solve this issue, we derive a novel fractional viscoelastic wave equation by combining an accurate relationship between angular frequency and complex wavenumber. The dissipation- and dispersion-dominated wave equations are also derived to simulate the amplitude-dissipation and phase-dispersion characteristics. The truncated Taylor-series expansion (TE) algorithm is developed to approximate the mixed-domain operators. After that, the generalized pseudospectral approach can be directly used to solve the new wave equation. In addition, an accurate viscoelastic wave equation constructed by the fractional time derivatives is used to calculate reference solutions to evaluate the accuracy of the new expression. Modelling results indicate that the newly proposed viscoelastic wave equation using the new fractional Laplacians is more accurate than the conventional one, especially in a small Q medium (i.e. QP = QS = 5). Furthermore, we also examine the accuracy of the TE approximation with a series of Q values. A homogeneous model and the modified BP2004 viscoelastic model are used to investigate the accuracy of viscoelastic wave propagations using the TE algorithm. All modelling results fully demonstrate the performance of the newly proposed viscoelastic wave equation and numerical algorithm. [ABSTRACT FROM AUTHOR]

Locating micro-earthquakes with high resolution and accuracy is a challenge for traveltime inversion, which has uncertainty on the order of a Fresnel zone (many wavelengths). We develop a wave-equation imaging method to increase resolution and reduce location errors to less than a wavelength, but requires very densely deployed receiver arrays with wide aperture and considerable computational cost. Instead of using acoustic data or direct P wave arrivals only, we use elastic multicomponent data and present a new method that uses the full P and S adjoint wavefields to image the microseismic source locations. We separate the P and S waves from the data, and extrapolate the P and S wavefields of each receiver subarray by solving the P and S adjoint wave equations in parallel. We formulate three source imaging conditions by multiplying over subarrays the adjoint P wavefield (I P), S wavefield (I S) and cross-correlated P and S wavefields (I PS). We perform numerical experiments on the highly realistic SEG SEAM4D reservoir model using surface acquisition array geometries. Results for 2-D and 3-D microseismic source estimations show clean images without noisy artefacts at shallow depths. In particular, I PS provides the highest resolution source location image, while I P is limited by the P wavelength and I S is influenced by small coda artefacts. The major-axis alignment and resolution of the source location image are determined by the hypocentral location with respect to the receiver array and illumination-angle coverage, respectively. We discuss the impacts of S -wave attenuation and frequency bandwidth on the source location images. Noise tests indicate that the imaging results are relatively insensitive to ambient noise, as is observed for the surface monitoring data. Using smoothed velocity models, the imaging results are similar to the results using the true realistically heterogeneous velocity model. The 90 per cent confidence ellipse of the source location due to Gaussian-distributed velocity errors shows a larger depth error as the source becomes deeper, while the horizontal error does not change as much. [ABSTRACT FROM AUTHOR]

To complement velocity distributions, seismic attenuation provides additional important information on fluid properties of hydrocarbon reservoirs in exploration seismology, as well as temperature distributions, partial melting, and water content within the crust and mantle in earthquake seismology. Full waveform inversion (FWI), as one of the state‐of‐the‐art seismic imaging techniques, can produce high‐resolution constraints for subsurface (an)elastic parameters by minimizing the difference between observed and predicted seismograms. Traditional waveform inversion for attenuation is commonly based on the standard‐linear‐solid (SLS) wave equation, in which case the quality factor (Q) has to be converted to stress and strain relaxation times. When using multiple attenuation mechanisms in the SLS method, it is difficult to directly estimate these relaxation time parameters. Based on a time domain complex‐valued viscoacoustic wave equation, we present an FWI framework for simultaneously estimating subsurface P wave velocity and attenuation distributions. Because Q is explicitly incorporated into the viscoacoustic wave equation, we directly derive P wave velocity and Q sensitivity kernels using the adjoint‐state method and simultaneously estimate their subsurface distributions. By analyzing the Gauss‐Newton Hessian, we observe strong interparameter crosstalk, especially the leakage from velocity to Q. We approximate the Hessian inverse using a preconditioned L‐BFGS method in viscoacoustic FWI, which enables us to successfully reduce interparameter crosstalk and produce accurate velocity and attenuation models. Numerical examples demonstrate the feasibility and robustness of the proposed method for simultaneously mapping complex velocity and Q distributions in the subsurface. Key Points: Subsurface velocity and attenuation structures are important for studying hydrocarbon reservoirs, geothermal dynamics, and magma chambersAn inversion scheme is developed for simultaneously mapping P wave velocity and attenuation based on a complex‐valued wave equationAn approximated Gauss‐Newton Hessian is utilized to reduce interparameter crosstalks and produce accurate velocity and attenuation models [ABSTRACT FROM AUTHOR]

Separating compressional and shear wavefields is an important step in elastic reverse-time migration, which can remove wave-mode crosstalk artefacts and improve imaging quality. In vertical (VTI) and titled (TTI) transversely isotropic media, the state-of-the-art techniques for wavefield separation are based on either non-stationary filter or low-rank approximation. They both require intensive Fourier transforms for models with strong heterogeneity. Based on the eigenform analysis, we develop an efficient pseudo-Helmholtz decomposition method for the VTI and TTI media, which produces vector P and S wavefields with the same amplitudes, phases and physical units as the input elastic wavefields. Starting from the elastic VTI wave equations, we first derive the analytical eigenvalues and eigenvectors, then use the Taylor expansion to approximate the square-root term in the eigenvalues, and finally obtain a zero-order and a first-order pseudo-Helmholtz decomposition operator. Because the zero-order operator is the true solution for the case of ϵ = δ, it produces accurate wavefield separation results for elliptical anisotropic media. The first-order separation operator is more accurate for non-elliptical anisotropy. Since the proposed pseudo-Helmholtz decomposition requires solving an anisotropic Poisson's equation, we propose two fast numerical solvers. One is based on the sparse lower-upper (LU) factorization, which can be repeatedly applied to the input elastic wavefields once computing the lower and upper triangular matrices. The second solver assumes the model parameters are laterally homogeneous within a given migration aperture. This assumption allows us to efficiently solve the anisotropic Poisson's equation in the z − kx domain, where kx and z denote the horizontal wavenumber and depth, respectively. Using the coordinate transform, we extend the pseudo-Helmholtz decomposition to the TTI media. The separated vector wavefields are used to produce PP and PS images by applying a dot-product imaging condition. Several numerical examples demonstrate the feasibility and applicability of the proposed methods. [ABSTRACT FROM AUTHOR]

In this paper, we derive a time-domain complex-valued wave equation for modelling visco-acoustic wave propagation in frequency-independent Q media. Starting from the frequency-domain visco-acoustic wave equation, we use a second-order polynomial to approximate the dispersion term and a pseudo-differential operator to approximate the dissipation term. These two approximations enable us to transform the frequency-domain visco-acoustic wave equation to the time domain. Due to the introduction of an imaginary unit in the dispersion approximation, the new wave equation is complex valued, which is similar to the time-dependent Schrödinger equation. The advantages of the proposed visco-acoustic wave equation include: (1) the dispersion and dissipation effects are naturally separated, which can be used to compensate amplitude loss in reverse-time migration by simply flipping the sign of the dissipation term; (2) the quality factor Q is explicitly incorporated in the wave equation, making it easier to derive the misfit gradient with respect to Q in full-waveform inversion in comparison with the traditional standard linear solid method; and (3) the new visco-acoustic wave equation can be numerically solved using finite-difference time marching and Fourier transform, which requires lower memory costs than frequency-domain visco-acoustic modelling solvers. The accuracy of the proposed equation is first verified by comparing with the analytical Green's function in a homogeneous medium. Then, two numerical examples are used to demonstrate that the new wave equation is capable of describing attenuation effects in heterogeneous frequency-independent Q media. [ABSTRACT FROM AUTHOR]

The acoustic behavior in fluid attenuating media can be effectively simulated using a fractional Zener model (FZM). Because of the fractional time derivatives of both stress and strain in the constitutive relationship, this mechanism is very realistic and flexible in describing seismic attenuation. However, using conventional FZM wave equations to propagate seismic waves requires storing large amounts of previous wavefield information to calculate the fractional time derivatives, which is unacceptable in practice. In this paper, we derive a new time-domain viscoacoustic wave equation in the framework of the FZM. This new equation does not contain any fractional time derivatives; thus, it is more economical in computational costs. Furthermore, the amplitude attenuation and phase dispersion effects are separated in the newly proposed equation, which is very favorable to compensate for energy loss and correct phase dispersion in reverse-time migration. To improve the accuracy, we incorporate a wave number (k)-space operator into the decoupled FZM wave equation to compensate for temporal dispersion errors caused by the second-order finite-difference discretization. Therefore, a high-temporal-accuracy viscoacoustic wave equation is derived to simulate nearly constant-Q wavefields in attenuating media. In the implementation, a low-rank decomposition method is introduced to solve the mixed-domain operators. Numerical analysis and modeling results demonstrate the effectiveness and applicability of the proposed method for simulating the decoupled viscoacoustic wavefield with high accuracy. [ABSTRACT FROM AUTHOR]

Compared with traditional adjoint-based migration, least-squares migration (LSM) can reduce finite-frequency effects, remove acquisition footprints and improve spatial resolution by solving a linear inverse problem for subsurface reflectivity. One important requirement for the success of LSM is having an accurate migration velocity model. Because of low signal-to-noise ratio (SNR), inaccurate traveltime picking, lack of low-frequency signals and limited acquisition aperture, it is still challenging to build an accurate velocity model using ray-based tomography or full waveform inversion. LSM with large velocity errors results in erroneous reflector locations, strong swing artifact and even non-convergence. To mitigate these issues, we develop a novel least-squares imaging framework in the subsurface half-opening angle domain. Instead of using high-wavenumber velocity perturbations as the reflectivity model as in traditional LSM, we parameterize the wave equation with an angle-dependent reflectivity, and derive the corresponding linearized forward modeling and adjoint migration operators. Because Gaussian Beam migration naturally incorporates propagation directions in wavefield extrapolation, we compute the Green's function using the Gaussian beam summation method. To improve the common-image gather (CIG) quality for low-fold and low-SNR data, a shaping regularization over the half-opening angles is introduced in the conjugate gradient scheme to iteratively update the angle-dependent reflectivity model. A flattening-enhanced summation is used to compute the stacked images by accounting for the depth moveout of CIGs caused by velocity errors, and produces constructive stacking results. Numerical experiments for benchmark models and a land survey demonstrate that the proposed method can improve LSM convergence and produce high-quality angle-dependent and stacked images even with inaccurate migration velocity models. [ABSTRACT FROM AUTHOR]