Research on seismic wave propagation in fluid-bearing reservoirs is crucial for unconventional oil and gas exploration. While traditional wave equation-based simulations are well-established, they often struggle to accurately model fluid-solid two-phase media characterized by complex tectonic features, strong discontinuities, and intricate interfacial couplings. This paper presents an advanced approach utilizing two mesoscopic discrete models: the lattice spring model (LSM) and the lattice Boltzmann model (LBM). The LSM simulates solid wavefields via the velocity Verlet algorithm, with optimized spatial distributions of springs and masses. Concurrently, the LBM is employed for fluid wavefield simulation, featuring a precisely designed lattice model and well-defined collision and streaming rules. In particular, a novel momentum exchange strategy incorporating the inverse squared distance method (ISDM) is constructed to enhance wavefield information transfer between the LSM and LBM at fluid-solid boundaries, aiming to improve accuracy and stability. The efficacy of the novel ISDM-LSM-LBM coupling algorithm is validated through simulations of layered models, fluid-filled cavity model, and river filling models. Compared to the original LSM-LBM algorithm, it exhibits enhanced computational accuracy and stability, albeit with a modest increase in computational burden and memory consumption. Comparative analysis with finite difference method (FDM) results reveals that the ISDM-LSM-LBM algorithm accurately captures wavefield characteristics. Although slight differences are observed in the depiction of reflected, transmitted, and converted waves at internal boundaries—attributable to the distinct fluid-solid interface treatments of each method—the new algorithm’s flexibility in boundary coupling and consideration of additional interaction forces offer significant improvements in accuracy and stability. This innovative simulation technique, independent of traditional wave equations, significantly enhances the precision and adaptability of wavefield calculations in fluid-solid media, providing substantial benefits for oil and gas geophysical exploration.

Full waveform inversion (FWI), recognized for its high precision, has garnered widespread application in velocity modeling. However, the computational cost of FWI continues to be a significant hurdle in its practical application. The FWI method employing multi-source simultaneous excitation, based on source-encoding, is a pivotal strategy for enhancing efficiency. However, conventional encoding methods suffer from crosstalk noise and are not adaptable to geometry where the positions of receivers are variable. To augment efficiency and ensure inversion accuracy, this paper introduces a time-frequency hybrid domain FWI algorithm that is founded on the resampled time-domain phase-sensitive detection (TD-PSD) method. This technique leverages the orthogonality of trigonometric integrals over complete periods to effectuate the decoupling of blended wavefields in the frequency domain and enables the computation of cross-correlation gradients without crosstalk. Concurrently, we propose a lossless sampling algorithm that markedly improves the efficiency of traditional TD-PSD integration calculations. Ultimately, the reliability of our proposed method is substantiated by numerical simulations on the two-dimensional (2D) overthrust model and three-dimensional (3D) field data, demonstrating that the inversion efficiency of our method is 4-6 times higher than that of traditional TD-PSD approach.

When using neural networks to predict low-frequency data of land seismic exploration, the biggest challenge is the difficulty in obtaining the wavelets of real seismic data. Training data that replicates the waveform structure of real seismic data are not easy to be generated. When a network trained with seismic records generated with other wavelets is directly applied to real data, it often can only approximate the prediction or, in some cases, cannot predict low-frequency data of land seismic data at all. Therefore, we propose a new approach that utilizes cross-convolutions to approximately unify the wavelets of real and simulated training seismic data. This allows the trained network to be used for predicting low-frequency data. Furthermore, the wavelets of the predicted low-frequency data are also known, which is more advantageous for subsequent full-waveform inversion. We employ two different wavelets to simulate real-world application scenarios. A simple layered model is utilized to validate the feasibility of the proposed method, and its generalization is tested on more complex Marmousi and overthrust models. For field data, better results can also be obtained using our methods.

Deconvolution technology, as an effective means to enhance the resolution of seismic data, has emerged as a prominent research area in the field of seismic exploration. However, due to its inherent ill-posed nature, seismic deconvolution poses significant challenges. The conventional approach for deconvolution employs a sparse inversion strategy to reconstruct underground reflection coefficients but overlooks the spatial relationship between adjacent seismic traces, resulting in inadequate spatial continuity of the deconvolved results. In this letter, we propose an adaptive seismic deconvolution strategy that incorporates spatial continuity regularization based on local seismic similarity. The adaptive multitrace deconvolution process consists of three parts: first, we consider spatial continuity by introducing a spatial regularization term derived from local seismic similarity; second, we formulate an objective function by combining$L1$-2 norm for sparse regularization with terms accounting for misfit and spatial regularization to reconstruct reflectivity; finally, we solve the objective function using alternative direction method of multipliers (ADMMs) and L-BFGS algorithms. Our proposed method effectively preserves weak effective signals while providing a clearer depiction of geological body distribution and ensuring superior spatial continuity in complex geological structures. Synthetic and field data tests demonstrate that our proposed method yields high-resolution deconvolved results with strong spatial continuity.

Joint estimation of elastic and petrophysical properties from seismic data and quantification of their uncertainties are critical aspects of reservoir characterization studies. The ensemble smoother with multiple data assimilation (ES-MDA) is proving to be a valuable tool for generating a reliable set of reservoir properties by matching simulated seismic responses with available observations. However, in standard implementations of ES-MDA, when the ensemble size is small, spurious correlations can arise in the cross covariances of model parameters and seismic data, leading to erroneous parameter updates in inappropriate regions. To mitigate this problem, this article introduces ES-MDA in conjunction with covariance localization (CL), termed ES-MDA_CL, which aims to reduce the impact of spurious covariance resulting from small ensembles and provide inversion results with robust error estimation. In the ES-MDA_CL framework, the model parameters are initially generated perturbatively by geostatistical simulations and subsequently updated while constrained by seismic data. The introduction of CL reduces the size of the initial ensembles, thereby increasing the computational efficiency of the entire inversion process. Through synthetic and field data tests, ES-MDA_CL demonstrates satisfactory inversion results with more reasonable computational times compared to ES-MDA. The proposed methodology enables the generation of inversion results with robust uncertainty estimation and holds promise for application to a wide range of geophysical challenges.

The complex-frequency-shifted convolutional perfectly matched layer (CFS-CPML) has been widely used in numerical simulations of the first-order seismic wave equations to avoid artifical reflections caused by truncated boundaries. However, numerically solving the second-order seismic wave equation is preferred for some large-scale geophysical algorithms, such as reverse-time migration (RTM). Extension of CFS-CPML to the second-order seismic wave equations has been realized by using a strict derivation approach in many literatures. However, the derived CFS-CPML formulations are complex, due to introduction of a large number of intermediate variables, which decreases overall computational efficiency, visibly. We present an efficient strategy to implement CFS-CPML in finite-difference time-domain (FDTD) modeling of the second-order seismic wave equations. We do not derive any continuous CFS-CPML formulations. Instead, we split the high-order centered-grid finite-difference (CGFD) stencil applied to approximate the second-order derivatives into two-level CGFD stencils. The inner stencils are viewed as the second-order CGFD approximations of the first-order derivatives at symmetric grid points. With this viewpoint, the convolution term related to CFS-CPML can be introduced to each inner CGFD approximation, naturally. The outer CGFD stencil is also viewed as CGFD approximation of a first-order derivative and it is augmented by a convolution term as well. By this way, CFS-CPML is incorporated into the FDTD simulations, efficiently. We take three-dimensional (3D) acoustic, viscoacoustic and elastic wavefields modeling for examples to verify the feasibility of the new implementation strategy of CFS-CPML.

The irregular surface topography has significant effects on seismic wave propagation, including introduction of free surface related multiples, strong scattering, and distortion of the shapes of seismic events. Because of the high efficiency, the finite-difference time-domain (FDTD) method is the most widely used numerical approach to solve three-dimensional (3D) seismic wave equations. However, FDTD at a uniform rectangular mesh could suffer from strong spurious diffractions, due to the stairgrid approximation to irregular surfaces. To resolve this problem, we immerse the irregular surfaces into fractional points under a uniform rectangular mesh and impose the free surface boundary condition along the normal direction of the surfaces. The implementation process of the free surface boundary condition in presence of irregular surfaces can be viewed a generalization of the traditional mirror image method designed for a planar free surface. The developed finite-difference (FD) scheme with the immersed boundary method (IBM), denoted as IBM-FD, needs to update the wavefields at a small number of ghost nodes, implicitly. We adopt Seidel iterations plus a third-degree Lagrange polynomial interpolation to realize this purpose. Compared with the FD scheme without IBM, the IBM-FD scheme only increases the computational cost slightly, due to the cost to determine the ghost nodes and update the wavefields at the ghost nodes. We present 3D scalar wavefield simulation examples with irregular surface topography to confirm the accuracy and stability of IBM-FD.

When the structure of the media is complex or there are strong interfaces, the simulated seismic waves obtained by the traditional forward modeling methods are often difficult to fit the real wavefields. To get more accurately simulated waves, in this article, we utilize the widely used 3-D lattice Boltzmann model (LBM) in fluid mechanics as a new tool for seismic wavefield simulation. In order to establish a general framework for the LBM in the 3-D case, first, we introduce an improved LBM scheme in 3-D form considering the reflection and transmission effects. Then, we carry out a systematic dispersion analysis of the LBM with different lattices by analogizing the dispersion derivation of the finite-difference method (FDM) and comparing the results of the LBM with those of the FDM. Finally, we explore the stability conditions and present the treatment of initial and boundary conditions. To verify the feasibility of this 3-D-LBM strategy, we perform 3-D wavefield simulation tests on a homogeneous model, a cavity model, and a five-layer model, respectively. By comparing the results obtained by LBM and FDM, we conclude that the accuracy of this improved LBM scheme is acceptable. This 3-D-LBM modeling scheme provides a new option for seismic wavefield simulation for complex media due to its flexibility of boundary processing and efficient parallelism capability.

Full-waveform inversion (FWI), as a high-precision model building method, plays an invaluable role in seismic exploration. The accuracy of conventional FWI is universally reduced by the cycle skipping, which can be improved by the optimal transport distance (OTD) objective function. However, the OTD-based FWI cannot converge to meaningful results with an inaccurately estimated source wavelet. To solve this dilemma, we construct a novel convolutional Wasserstein distance (CW) objective function by applying the OTD objective function to convolved seismograms. Before the standard non-negative and normalization preprocessing of OTD, we first convolve the observed data with a reference trace selected from simulated seismograms and convolve the simulated data with a trace selected from the observed data. The both convolved data sets are naturally regarded as with an identical source, so the data difference caused by the inaccurately estimated source wavelet is eliminated. The adjoint source corresponding to the new objective function is derived. The velocity model can be updated by using a quasi-Newton method according to the gradient of the objective function generated by the adjoint-state method. We investigate the effectiveness of our objective function by one-dimensional signals and several FWI examples. Furthermore, this new objective function still delivers an excellent performance in releasing the local minimum problem and resisting noise when the wavelet used in FWI is inaccurate.

Staircase approximation in the simulation of the arbitrarily irregular surface will cause artificial interference of wavefield. The body-fit coordinate transformation can accurately simulate the arbitrarily irregular surface and avoid the staircase approximation. The finite difference method (FDM) based on body-fit coordinate system (BFCS) can simulate the propagation process of wavefield in the media with an undulating surface. Currently, the collocated grid is usually used in BFCS. However, adopting this method to simulate waves based on the first-order velocity-stress wave equation in the BFCS will cause oscillation due to the odd-even decoupling; therefore, additional filtering operations are required. To solve the above problem, we propose collocate-staggered grid (CSG) to simulate the first-order velocity-stress wave equation in the 2-D BFCS. CSG uses staggered grid scheme to calculate the equation parameters and adopts collocated grid scheme to store these values. We use second-order accuracy spatial and temporal difference for numerical tests. The results show that the elastic wavefield modeling method in the BFCS using CSG difference scheme has the characteristics of high accuracy and stable. This method is suitable for the wavefield simulation in a medium with a large Poisson’s ratio.

Full-waveform inversion (FWI) serves as a useful tool to quantitatively investigate the properties of the subsurface. Presently, 3-D elastic FWI uses a finite-difference time-domain (FDTD) approach in numerical simulation. However, such an FDTD scheme often includes only second-order temporal approximations, causing errors in temporal dispersion in the case of a large time-stepping size. Such temporal dispersion will affect the inversion results and reduce the inversion quality. We introduce a unique 3-D elastic FWI using a temporal fourth-order finite-difference (FD) approximation. A new quasi-stress–velocity elastic equation is solved by the temporal fourth-order and spatial arbitrary even-order FDTD method, and a novel inversion procedure for the convolutional objective function based on this equation is derived. The multiscale strategy is used to enhance the robustness of our algorithm. The forward modeling and FWI examples presented here demonstrate that our method can achieve modeling and inversion with a high degree of accuracy.

Absorption in subsurface media severely degrades seismic data quality. Seismic attenuation compensation as an important processing method can effectively improve the resolution and fidelity of seismic data. Based on sparse reflectivity model and attenuated convolution function, inversion-based compensation approaches show better stability and accuracy over traditional direct compensation schemes. However, conventional inversion-based compensation methods are conducted on single trace, which ignore the subsurface spatial continuity and make the compensated result contaminated with high-frequency noise. In this paper, we develop a structurally constrained multichannel L1-2 minimization for seismic attenuation compensation. We first estimate structure tensors from migrated seismic images. The structure tensors can be decomposed by eigenvalues and eigenvectors, which can reflect the structural orientations. Then, we introduce the estimated orientations as a regularization term to the L1-2 inversion-based compensation objective function. In this way, we can improve the stability of the compensation result and enhance the spatial continuity of the compensated seismic reflectors. The structure-guided L1-2 regularized compensation objective function can be efficiently solved via difference of convex algorithm and alternating direction method of multipliers. Synthetic and field data examples demonstrate that the proposed method possesses superior performance over conventional L1-2 regularized inversion-based compensation.

Reverse-time migration (RTM) has the particular capacity for complex geological structure imaging. However, massive storage and high computational costs caused by conventional cross-correlation imaging conditions restrict the large-scale application of RTM. The excitation amplitude imaging condition (EAIC) has the cheapest imaging cost, but inadequate wavefield information causes less tolerance for noise. To alleviate these limitations, we introduce a local Nyquist cross-correlation imaging condition, which serves as a transition between the conventional cross-correlation imaging conditions and the EAIC. Instead of using the full wavefields or only the excitation amplitude to construct imaging, the local cross-correlation imaging condition (LCIC) utilizes the wavefields around the corresponding time of the maximum amplitude at each grid during the source wavefield simulation. Moreover, considering the possible oversampling situations in the LCIC, the Nyquist sampling theorem is flexibly embedded into the local cross-correlation scheme to further reduce the storage costs and improve the efficiency of RTM. Compared with the cross-correlation imaging conditions, the proposed strategy can obtain the similar results and reduce the storage requirement and time costs significantly. In the meantime, it maintains better adaptability to the complex imaging environments than the EAIC. Migration tests of synthetic and field datasets demonstrate that the local Nyquist cross-correlation scheme features good feasibility, efficiency, and practicability in RTM. As a consequence, the proposed local Nyquist cross-correlation imaging condition can effectively save storage and time costs and provide a reliable migrated image even from noisy observed data.

We develop a fractional Laplacian viscoacoustic wave equation-based full waveform inversion (FWI) method. The main novelty is efficient reconstruction of the source wavefields in gradient computation by the adjoint state method. Our FWI is based on Fourier pseudo-spectral time-domain (PSTD) numerical solutions of the fractional Laplacian viscoacoustic wave equation, which can describe the frequency independent${Q}$(quality factor) behaviors of seismic waves accurately. The presented wavefield reconstruction strategy utilizes reverse time-marching formulae to compute the source wavefields backward in time, including an implicit formula in the interior domain and an explicit one in the perfectly matched layer (PML) absorbing domains. Since the wavefields in the entire domain are recovered, our method does not require storing massive boundary values like conventional reconstruction methods. To avoid numerical instability in reverse time reconstruction, we design a global-local checkpointing technique. When reverse time reconstruction encounters numerical instability, the forward propagation restarts from the nearest global checkpoint and proceeds to the unstable time step. Then, the reverse time reconstruction continues. The local checkpoints help to delay the numerical instability in the PML domains. Numerical examples verify the feasibility of our reconstruction strategy and the efficiency gain achieved by using this strategy in FWI.

Wave simulation in absorptive media using decoupled fractional Laplacian wave equation has received widespread attention in recent years, largely due to its precise description of frequency independent $Q$ and easy attenuation compensation in seismic processing. With many algorithms to solve the fractional Laplacian, $k$ -space pseudo-spectral method is predominantly used in academia, where the computing facilities support fast Fourier transforms. However, its global nature prevents the parallelization and computational efficiency of the forward solver, especially for large-scale applications. We propose to solve viscoacoustic wave equation using domain decomposition. A local Fourier basis is constructed around the truncated area to improve the periodicity and smoothness of the decomposed wave information. After independently simulating in the subdomains, a pointwise patching procedure is applied to maintain the continuity of the wavefield between subdomains. Numerical experiments show that this new algorithm obtains high computational efficiency without compromising the numerical stability condition of the traditional pseudo-spectral method. This work becomes more attractive for seismic inversion and imaging problems by improving its parallelization.

Spikelike noise is a common type of random noise existing in many geoscience and remote sensing data sets. The attenuation of spike-like noise has become extremely important recently, because it is the main bottleneck when processing the simultaneous source data that are generated from the modern seismic acquisition. In this letter, we propose a novel low-rank decomposition algorithm that is effective in rejecting the spike-like noise in the seismic data set. The specialty of the low-rank decomposition algorithm is that it is applied along the morphological direction of the seismic data sets with a prior knowledge of the morphology of the seismic data, which we call local slope. The seismic data are of much lower rank along the morphological direction than along the space direction. The morphology of the seismic data (local slope) is obtained via a robust plane-wave destruction method. We use two simulated field data examples to illustrate the algorithm workflow and its effective performance.