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.

Bayesian inversion is capable of integrating seismic data, well-log data, and geological data to obtain a posterior probability distribution function (PPDF) of elastic parameters. Due to the absorption of strata, amplitude attenuation and phase distortion inevitably occur during the propagation of seismic wave, resulting in low resolution of seismic data. The traditional linear Bayesian inversion method is based on the stationary convolution model, and amplitude compensation and phase correction are needed to be conducted in advance. Uncertainties exist in both compensation and inversion, and the inversion results cannot account for the uncertainties in the directly observed seismic data resulting in accurate inversion results. To estimate uncertainty more accurately and improve inversion accuracy, in this article, we develop a nonstationary prestack linear Bayesian stochastic inversion (NSPLBSI) method. Through the proposed method, prior information from well-log data can be effectively introduced to constrain inversion. Compensation and inversion are integrated into one procedure, which can estimate the posterior uncertainties more accurately directly from seismic data compared with traditional two-step inversion methods, i.e., first compensating attenuation and second performing inversion. Also, more accurate inversion results are obtained. In order to verify the rationality and effectiveness of our proposed method, we conduct inversions using synthetic seismic data on the Marmousi model and a field seismic dataset. Numerical examples show that it can not only compensate amplitude well but also obtain high-precision inversion results with smaller prediction interval.

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.

High-resolution prestack inversion combining the well-logging and seismic data is a significant geophysical task and can be achieved by two kinds of stochastic inversion approaches, the geostatistical inversion (GSI) and Bayesian linearized inversion (BLI). Nevertheless, the existing GSI is restricted by the heavy iteration calculation. Although BLI can avoid this issue, it suffers from the large core matrix inverse. A fast multi-Gaussian inversion (FMG_INV) is proposed herein to achieve the well-log and seismic combined inversion with higher efficiency than GSI and BLI. FMG_INV is derived from prestack BLI, which requires a large core matrix inverse. However, FMG_INV utilizes a simplification strategy and reduces the core matrix dimension of BLI. This improvement is presented under the assumption of statistical independence between well-logging and seismic data, which relieves the issue of large matrix inverse in BLI to a great extent. Moreover, the spatial and statistical correlation between different parameters in prestack stochastic inversion is presented by a multi-Gaussian distribution and may reduce inversion accuracy, and FMG_INV solves this problem by a novel decorrelation strategy. The 1-D, 2-D, and 3-D field tests and a synthetic data test are given herein to verify the effectiveness of FMG_INV. The 1-D and 2-D tests of traditional BLI are also conducted for comparison. The results demonstrate that FMG_INV achieves the same satisfying inversion accuracy and resolution with BLI but much lower time consumption than BLI.

Bayesian linearized inversion (BLI) stands out as an exceptional stochastic inversion method in the realms of geophysics and remote sensing. It excels in estimating inversion results and assessing their uncertainty with remarkable efficiency. However, one of the challenges faced by BLI lies in the inversion of its core matrix. To surmount this limitation, an innovative dimension reduction strategy is proposed based on the discrete cosine transform (DCT), thus formulating a rapid BLI approach termed DCT-BLI. Within this method, the DCT-based reduction strategy effectively compresses a large sparse matrix by extracting its essential information, transforming the inversion of this sizable matrix into the inversion of a reduced-size counterpart. A compression factor (CF), defined as the size ratio of matrices after and before reduction, quantifies the extent of matrix reduction. DCT-BLI integrates the strengths of both BLI and the DCT-based reduction strategy. Leveraging this reduction approach, DCT-BLI tackles the challenge of inverting its sizable core matrix. Through the synthetic and field data tests, DCT-BLI exhibits clear superiority over BLI in terms of efficiency, and the DCT-based reduction method achieves a remarkable two-thirds reduction in the core matrix size of BLI without compromising inversion accuracy.

Time-domain elastic full waveform inversion (FWI) uses seismic data to recover the high-resolution subsurface medium properties for structural imaging and lithologic identification. The approximate pulse generated by the explosion source in seismic exploration evolves into a limited bandwidth seismic wavelet through the propagation of the seismic wave in the underground medium, exhibiting strong energy near the dominant frequency and weak energy far away. Therefore, time-domain FWI using the band-limited seismic wavelet can only match the energy of the dominant frequencies present in a dataset to a large extent, resulting in insufficient low-wavenumber updates of model parameters because of the weak energy of low frequencies. To remove the effect of finite-frequency-band wavelet spectra from the time-domain FWI, we propose a time-domain elastic FWI without wavelet spectral limitation. In our method, a newly refined seismic wavelet with a normalized amplitude and accurate phase is used to propagate seismic waves in the time domain. Data residual measurement is performed based on the summation of single-frequency residuals between the normalized synthetic and observed frequencies. The adjoint-state method is used to approximate the gradients of the model parameters, and the decoupled wavefields obtained by the phase-sensitive detection method were involved in the gradient calculation. Overall, the proposed method helps FWI to avoid falling into local minima by enhancing the low-wavenumber reconstruction of the velocity models. The elastic FWI numerical tests and field data FWI application demonstrate that our approach can reliably recover high-precision inversion results.

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.

Lattice spring model (LSM) is a novel method to simulate seismic wave propagation from a micromechanical perspective, which describes the elastic dynamics in complex media comprehensively and delineates the dynamic characteristics of the wavefield. However, it remains a great challenge to simulate wavefields at irregular interfaces with high accuracy. This work presents an adaptive location LSM (ALLSM) method to substantially improve the accuracy of simulated wavefield in a solid model. First, we design an algorithm to pick up the intersection coordinates of the internal interfaces with regular mesh. We obtain the intersection coordinates by solving the defined interface functions with different kinds of spring functions in the LSM. Then, we calculate the weights of springs with different elastic coefficients in the whole grid by the intersection coordinates. Finally, we obtain the real spring elastic coefficients in the presence of irregular interfaces by the arithmetic weighted average method, and develop the corresponding interface position ALLSM to simulate the wavefield propagation mechanism, in order to substantially improve the accuracy of simulation. Numerical tests show that the improved ALLSM is more accurate than conventional LSM and second-order finite difference method (FDM). The accuracy of simulated wavefields is significantly improved about 20%. The proposed method offers a new tool for wavefields simulation in complex media.

The establishment of initial subsurface model is a crucial step for seismic inversion. An accurate and reasonable initial model can mitigate the ill-posedness of seismic inversion and improve the quality of inversion result. A common method for building initial model is the well-log data interpolation. However, the traditional well-log interpolation method ignores the structural information of the subsurface, resulting in the constructed initial model lacking geological meaning. We propose a novel structurally constrained modeling method (SCMM) to obtain a geologically reasonable initial impedance model for poststack seismic inversion. Well-log interpolation can be represented as an inverse problem. SCMM constrains the inversion process by using a regularization operator that forces the well-log data to be extended to the entire seismic working area along the subsurface local structural direction. First, we calculate the seismic dip from the poststack seismic profile. Then, we design the structural operator based on the estimated seismic dip information to constrain the interpolation process. Under the framework of inversion, the interpolation objective function can be established by combining the structural operator with the well-log data misfit term, and it can be solved efficiently by the conjugate gradient algorithm. Synthetic and field data tests show that the initial model built by SCMM is more consistent with the geological rules than that built by traditional method, and the poststack impedance inversion using SCMM is better than that using traditional modeling method in terms of convergence property and accuracy of inversion result.

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.

Impedance inversion is an effective method to estimate the properties of subsurface model from poststack seismic data. However, owing to the low quality of seismic data and other reasons, impedance inversion methods usually suffer from spatial discontinuities and low resolution. In order to overcome these problems, we have developed an impedance inversion method based on geological structure-guided total variation (GSGTV) constraint. Different from traditional total variation (TV) constraint, GSGTV considers the spatial distribution of geological structures instead of imposing gradient constraints on impedance only in fixed directions (horizontal and vertical). Therefore, GSGTV not only retains the advantage of traditional TV that can preserve the edge of layer but also overcomes the shortcomings of traditional TV that is only suitable for inverting block structures. This ensures that GSGTV can simultaneously improve the resolution and spatial continuity of the inversion result. The realization of impedance inversion method based on GSGTV constraint requires three steps. First, estimate the local structural orientation of the subsurface from seismic data. Then, construct GSGTV regularization term based on the local structural orientation. Finally, establish the inversion objective function and solve the function by the split-Bregman iterative algorithm. We compared the proposed method and the impedance inversion method based on traditional TV constraint with synthetic and real seismic data. The inversion results confirm that our method can improve the resolution and spatial continuity.

The increasing complexity of seismic exploration projects and the request for higher imaging resolution have driven the geophysics community to look for a sound understanding of the subsurface formation to optimize seismic structure interpretation and reservoir characterization. Crosswell seismic survey aims at obtaining higher resolution images of the interwell regions and more accurately characterizing the reservoir dynamics. However, the presence of the intrinsic seismic attenuation of rocks as seismic waves propagate through the subsurface results in amplitude decay and velocity dispersion. This inevitably decreases the imaging resolution and the reliability of the subsequent seismic interpretation and reservoir characterization. To compensate for the attenuating effect, one may restore to attenuation compensation technique during seismic imaging. We here present the$Q$-compensated viscoelastic reverse time migration ($Q$-ERTM) based on the decoupled fractional Laplacian (DFL) viscoelastic wave equation for high-resolution crosswell imaging. We develop an explicit stabilization scheme to resolve the cumbersome numerical instability issue in$Q$-ERTM. The merits of explicit stabilization are twofold. First, it simplifies the workflows of the$Q$-ERTM by avoiding domain transforms. In addition, it provides a flexible way for stabilization parameter tuning by introducing a reference scaling factor. We follow the best practices of high-performance computing with the MPI + CUDA configuration for numerical implementation. A toy crosswell imaging example and a more realistic time-lapse crosswell seismic survey with a CO2plume injection are provided to verify the feasibility and stability of the proposed method.

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.

Seismic deconvolution is an effective approach to improve the resolution of seismic data and plays an important role in migration imaging, reservoir prediction and other fields. However, conventional deconvolution methods are usually based on sparse-type regularization (e.g.,$L_{1}$-norm) and adopt a trace-by-trace inversion strategy to reconstruct the subsurface reflectivity series. Although such methods can improve the resolution of seismic records to a certain extent, the lack of spatial constraint will result in poor spatial continuity in the reconstructed reflectivity. This phenomenon is particularly obvious in regions with complicated geologic structures. For the purpose of overcoming this issue, we have developed a structure-oriented regularization-based multichannel sparse spike deconvolution (SOR-based MSSD) method. This method imposes$L_{1}$-norm regularization on the reflectivity to obtain the high-resolution subsurface reflectivity series and imposes structure-oriented regularization (SOR) on the expected high-resolution seismic data to improve the spatial continuity of the inversion result. First, we construct SOR term based on the local structural orientations which can be estimated from the poststack seismic data. Then, we integrate the seismic data misfit term, the$L_{1}$-norm constraint term, and the SOR term to formulate the objective function. At last, we use the alternating direction method of multipliers (ADMMs) to efficiently solve the objective function. We compare the SOR-based MSSD with existing methods by using synthetic and field data. Both deconvolution examples illustrate the performance of proposed method in terms of improving the spatial continuity.

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.

Seismic deconvolution is a typical ill-posed inverse problem. The regularization technique in terms of different prior information is used for a unique and stable solution. Due to the difference between prior information and the actual subsurface situation, it is hard to obtain a solution with satisfactory accuracy and resolution. This letter presents the dictionary learning as an efficient adaptive deconvolution method for the reflectivity reconstruction problem. Considering the curse of dimensionality of conventional dictionary learning and the suboptimal solution of the patch-based dictionary learning, we take the convolutional sparse coding (CSC) model as the dictionary learning method. In this method, the prior information can be obtained from the well-log data in the form of sparse CSC dictionary of reflectivity. On the assumption that the deposition of the subsurface layers is stable, the CSC dictionary extracted from the well-log data can also be applied in the whole work area. The CSC-based deconvolution can be seen as the adaptive deconvolution due to the independence of the assumption made about the reflectivity and seismic data. The process of the adaptive CSC-based deconvolution is divided into three parts. First, the CSC dictionary is learned from the well-log data. Then, the objective function is formulated by combining the CSC dictionary and the single-channel seismic record misfit term for the reconstruction of reflectivity. Finally, the objective function is efficiently solved with the coordinate descent approach. We illustrate the performance of our adaptive deconvolution with synthetic and field seismic data.

Simultaneous-source acquisition, breaking the limit of conventional seismic acquisition, is a rapidly evolving research field, due to its advantage in reducing survey time and improving data quality. The benefits of simultaneous-source acquisition are compromised by the intense blending interference. Separating a blended record into a group of individual records, known as “deblending” is one of the most popular solution to the problem. However, the blended records are often corrupted by random noise, which causes difficulties in separation. In an iterative deblending algorithm, the incoherent interference can be simulated and subtracted from the blended record. When the random noise is strong, it is difficult to simulate the incoherent interference. In this paper, we propose a hybrid-sparsity constraint model that applies the dictionary learning into the deblending framework that is based on the sparsity-promoting transform to deal with extremely noisy simultaneous source data. The dictionary learning with fine-tuned adaptation can learn the incoherent interference into atoms and reject random noise. Then, the sparse transform-based framework is implemented to iteratively separate the signal and interference. We use two synthetic examples to demonstrate the advantage of the proposed method in extremely noisy situations. Two field examples further confirm the superior deblending performance of the proposed method for the noisy simultaneous-source data over the curvelet transform-based and rank reduction-based methods.

The proximal splitting algorithm, which reduces complex convex optimization problems into a series of smaller subproblems and spreads the projection operator onto a convex set into the proximity operator of a convex function, has recently been introduced in the area of signal processing. Following the splitting framework, we propose a novel three-operator proximal splitting (TOPS) algorithm for 3-D seismic data reconstruction with both singular value decomposition (SVD)-based low-rank constraint and curvelet-domain sparsity constraint. Compared with the well-known forward-backward splitting (FBS) method, our proposed TOPS algorithm can be flexibly employed to recover a signal satisfying double convex constraints simultaneously, such as low-rank constraint and sparsity constraint used in this letter. We have used both synthetic and field data examples to demonstrate the superior performance of the TOPS method over traditional SVD-based low-rank method and curvelet-domain sparsity method based on the FBS framework.

The traditional high-order staggered-grid finite-difference (SGFD) method has high-order accuracy in space, but only the second-order accuracy in time, which makes the traditional SGFD method suffer from a large temporal dispersion error during long-distance wave propagation. This paper develops temporal fourth- and sixth-order and spatial arbitrary evenorder SGFD schemes to model isotropic elastic wave propagation. The temporal high-order SGFD schemes have smaller temporal dispersion than the traditional temporal second-order scheme, and thus allow larger time steps to attain a similar accuracy. The developed temporal high-order SGFD schemes are applied to simulate a quasi-stress–velocity wave equation (QWE) that is derived in the framework of a$k$-space approach. A split QWE (SQWE) is further developed, and numerical simulation of SQWE results in separated P (compressional)-wave and S (shear)-wave. Theoretical computational cost analysis verifies that the numerical simulation of QWE using the temporal fourthand sixth-order SGFD schemes is more efficient than the numerical simulation of the traditional stress–velocity wave equation using the traditional temporal second-order SGFD scheme in 2-D. In 3-D, the temporal fourth-order SGFD scheme still runs faster than the traditional temporal second-order scheme; however, the temporal sixth-order scheme is more efficient only when a longer stencil length than 12 is adopted. Numerical examples confirm the correctness of the developed elastic wave modeling schemes.

Seismic data interpolation or reconstruction plays an important role in seismic data processing. Many processing steps, such as high resolution processing, wave-equation migration, amplitude-versus-offset and amplitude-versus-azimuth analysis, require regularly sampled data. The reconstruction can be posed as an inverse problem, which is known to be ill posed and requires constraints to obtain unique and stable solutions. In this letter, we propose an iterative scheme to interpolate the big gaps with a slope constraint. In the first iteration, the smooth radius must be large to estimate the smooth dip from the decimated data, and a large scaling parameter can guarantee the stability of the inversion. In the later iterations, the smooth radius will be shortened in order to get a more accurate dip estimation from the updated result. When the dip estimation is accurate, a small scaling parameter can not only guarantee the convergence of the inversion but also obtain a result with high signal-to-noise ratio. We compare the proposed method with the well-known projection-onto-convex-sets method on synthetic and field data examples. The interpolation results illustrate the advantage of the proposed method in interpolating the big gaps.