Search

Full-waveform inversion (FWI) is an effective method that uses kinematic and dynamic information to provide an accurate image of the subsurface target. By applying the Bayesian approaches to FWI workflow, the traditional deterministic solution is replaced with a maximum a posteriori solution. This provides not only a reliable subsurface model for the nonlinear inverse scattering problem but also an uncertainty quantification. The high computational cost of the Bayesian FWI and the complexity of a large-scale anisotropy inversion, however, restrict its application in industrial settings. We have developed an FWI approach for tilted transversely isotropic (TTI) media. The strain field and the particle displacement of the anisotropic media are expressed by an integral method based on the scattering theory and an explicit multiscattering modified Green’s function. To reduce the computational cost, the wavefield is calculated by using the preconditioned generalized minimal residual iterative solver, which approximates the solution in the Krylov subspace. The Fréchet derivatives are expressed by using the distorted Born iterative method. An iterated extended Kalman filter is applied to solve the inverse problem and the posterior covariance matrix frequency-by-frequency. Numerical experiments indicate that the presented methodology can reasonably reconstruct the elastic moduli of the TTI medium. Furthermore, the numerical experiments demonstrate the accuracy and feasibility of our methodology.

Correct subsalt structure recovery is important for subsalt reservoirs exploration. However, high precision subsalt structure retrieving using full waveform inversion (FWI) is not easy caused by the shielding from salt body and the depth of the target. In order to get correct subsalt structures, techniques which can effectively work with deep targets are desired. Besides, correct recovery of the salt body is critical for the recovery of subsalt structure. However, the salt body recovery itself is also challenging in FWI, especially when a reasonable starting model is not available. In this work, we consider both the recovery of the salt structure and subsalt structure. A local angle domain illumination compensated direct envelope inversion with local resolution function is proposed to provide large-scale background salt structure with enhanced and more balanced recovery of the salt body in deep part. Then FWI with angle domain illumination compensation is used to improve the recovery for both the salt and subsalt structures. To further enhance the quality of the inversion, model updating with L-BFGS is embed into the whole process. Numerical examples with the SEG/EAGE salt model and the Sigsbee model demonstrated that the proposed method can retrieve the salt structure as well as the subsalt structure effectively with high precision starting from a linear starting model.

Joint inversion algorithms continue to receive a lot of interest in the application of geophysical interpretation due to the reduction of uncertainty of the inverse solutions. Over the past decade, several coupling constraints, structural or rock physics, have been introduced to the joint inversion of multiple geophysical data sets. However, a single coupling constraint may not provide the optimum solution for all geologic scenarios. To improve the application of the multiple geophysical data sets, in which some rock-physics information is available, we have developed a new joint inversion algorithm based on a mixed structural and rock-physics coupling constraints. The presented algorithm incorporates the cosine dot product gradient function, a structural coupling constraint, with a guided fuzzy c-means clustering strategy, which is a rock-physics coupling constraint. The former imposes the structural similarity among different model parameters, whereas the latter improves the linear or nonlinear correlation among model parameters. We use the exact structural resemblance algorithm to numerically solve the mixed coupling constraint joint inversion objective function. The developed algorithm is applied on 2D synthetic and field data examples. Our results find a significant improvement in the resolution of the reconstructed models, in comparison with traditional single coupling algorithms, even when the rock-physics information is not complete. This algorithm is more conducive to obtaining clear geologic unit boundaries and provides a favorable basis for delineating metallogenic target areas.

In seismic data processing, compensation for geometric spreading along the raypath is necessary. The accurate approximation of the relative geometric spreading is critical because the traveltime parameters in the approximation can be estimated from velocity analysis. Two approximation types (direct and indirect) are generalized based on the procedure to derive the explicit expression. We develop perturbation-based approximations of the relative geometric spreading for P and S waves with these two types in the elastic transversely isotropic media with a vertical symmetry axis (VTI). The rational approximation and the general moveout approximation counterparts are compared to verify the superiority of our approach. Our study compares the performance of direct and indirect approximations for optimal selection under different cases. In the numerical examples, we apply the error sensitivity in the layer depth, anellipticity value, and vertical velocity ratio. It is found that the P-wave case is much more accurate than the S-wave case, and the direct type obtains relatively higher accuracy than the indirect one. In addition, the perturbation coefficients for the acoustic VTI model are derived not only with a much simpler expression but also achieve a highly accurate result that is recommended for the P-wave seismic processing.

Large-scale geologic structures with strong contrasts present difficulties for seismic imaging and inversion. Previously, direct envelope inversion (DEI) with ultralow frequency has been proposed to retrieve strong-contrast velocity models. However, the envelope signals lose the instantaneous phase of the original wavefield, resulting in erroneous seismic inversion results. To solve this problem, the instantaneous phase of the wavefield is encoded into envelopes for DEI. In this way, a polarized envelope is created to contain both abundant low frequencies and the instantaneous phase information of the original seismic waveform. Furthermore, to improve the resolution of deep regions underneath the strong reflections, we propose a phase-amplitude-based polarized direct envelope inversion (PA-PDEI) misfit in the time-frequency domain that uses an amplitude factor to emphasize the phase information. This makes it possible to adjust the amplitude influence and boost the low-amplitude seismic signals for deep targets. Tests on a model with a salt body and on field data show that combining the PA-PDEI with the phase-amplitude-based full-waveform inversion (PA-FWI) make it possible to reliably reconstruct velocity models, especially for deep regions, despite large and strong perturbations.

Simulation of the seismoelectric effect serves as a useful tool to capture the observed seismoelectric conversion phenomenon in porous media, thus offering promising potential in underground exploration activities to detect pore fluids such as water, oil, and gas. The static electromagnetic (EM) approximation is among the most widely used methods for numerical simulation of the seismoelectric responses. However, the static approximation ignores the accompanying electric field generated by the shear wave, resulting in considerable errors when compared to analytical results, particularly under high-salinity conditions. To mitigate this problem, we have adopted a spatial high-order finite-difference time-domain method based on Maxwell’s full equations of time-varying EM fields to simulate the seismoelectric response in 2D mode. To improve the computational efficiency influenced by the velocity differences between seismic and EM waves, different time steps are set according to the stability conditions, and the seismic feedback values of EM time nodes are obtained by linear approximation within the seismic unit time step. To improve the simulation accuracy of the seismoelectric response with the time-varying EM calculation method, finite-difference coefficients are obtained by solving the spatial high-order difference approximation based on the Taylor expansion. Our method yields consistent simulation results compared to those obtained from the analytical method under different salinity conditions, thus indicating its validity for simulating seismoelectric responses in porous media. We further apply our method to layered and anomalous body models and extend our algorithm to three dimensions. Results indicate that the time-varying EM calculation method can effectively capture the reflection and transmission phenomena of the seismic and EM wavefields at the interfaces of contrasting media. This may allow for the identification of abnormal locations, thus highlighting the capability of seismoelectric response simulation to detect subsurface properties.

Many applications in seismology involve the modeling of seismic-wave traveltimes in anisotropic media. We have developed homotopy solutions of the acoustic eikonal equation for P-wave traveltimes in attenuating transversely isotropic media with a vertical symmetry axis. Instead of the commonly used perturbation theory, we use the homotopy analysis method to express the traveltimes by a Taylor series expansion over powers of an embedding parameter. For the derivation, we first perform homotopy analysis of the eikonal equation and derive the linearized ordinary differential equations for the coefficients of the Taylor series expansion. Then, we obtain the homotopy solutions for the traveltimes by solving the linearized ordinary differential equations. Results of our investigation with approximate formulas demonstrate that the analytical expressions are efficient methods for the computation of traveltimes from the eikonal equation. In addition, these formulas are also effective methods for benchmarking approximate numerical solutions in strongly attenuating anisotropic media.

The conventional coherence attribute is typically applied to migrated full-stack seismic data volumes to detect geologic discontinuities. Recently, multispectral, multiazimuth, and multioffset coherence attributes have been proposed and implemented with different seismic data volumes of specific frequencies, azimuths, and offsets to enhance discontinuities. Generally, geologic anomalies, such as faults and channels, will be better illuminated by a perpendicular rather than a parallel direction for computation. Therefore, we have developed a multidirectional eigenvalue-based coherence attribute by establishing multiple covariance matrices along certain different directions on a single poststack volume. We adopt two methods to compute a multidirectional coherence attribute. One is to compute multiple coherence volumes in different directions and to define the minimum as the final multidirectional coherence. This method is time-consuming, but it could provide partial and overall discontinuity simultaneously. The other method obtains one coherence volume by summing covariance matrices in different directions, which is computationally efficient, but only provides overall discontinuity. The performance of our 3D physical model and field data volumes demonstrates that multidirectional coherence can highlight subtle geologic structures with a higher resolution than conventional coherence. This suggests that a multidirectional coherence attribute may serve as an effective tool for detecting the distribution of geologic discontinuities in seismic interpretation.

We have developed a finite-difference iterative solver of the Helmholtz equation for seismic modeling and inversion in the frequency domain. The iterative solver involves the shifted Laplacian operator and two-level preconditioners. It is based on the application of the preconditioners to the Krylov subspace stabilized biconjugate gradient method. A critical factor for the iterative solver is the introduction of a new preconditioner into the Krylov subspace iteration method to solve the linear equation system resulting from the discretization of the Helmholtz equation. This new preconditioner is based on a reformulation of an integral equation-based convergent Born series for the Lippmann-Schwinger equation to an equivalent differential equation. We have determined that our iterative solver combined with the novel preconditioner when incorporated with the finite-difference method accelerates the convergence of the Krylov subspace iteration method for frequency-domain seismic wave modeling. A comparison of a direct solver, a one-level Krylov subspace iterative solver, and our two-level iterative solver verified the accuracy and accelerated convergence of the new scheme. Extensive tests in full-waveform inversion demonstrate the solver’s applicability to such problems.

As an approximation of the spherical-wave reflection coefficient (SRC), the plane-wave reflection coefficient does not fully describe the reflection phenomenon of a seismic wave generated by a point source. The applications of SRC to improve analyses of seismic data have also been studied. However, most of the studies focus on the time-domain SRC and its benefit to using the long-offset information instead of the dependency of SRC on frequency. Consequently, we have investigated and accounted for the frequency-dependent spherical-wave reflection coefficient (FSRC) and analyzed the feasibility of this type of inversion. Our inversion strategy requires a single incident angle using reflection data for inverting the density and velocity ratios, which is distinctly different from conventional inversion methods using amplitude variation with offset. Hence, this investigation provides an alternative approach for estimating media properties in some contexts, especially when the range of aperture of the reflection angles is limited. We apply the FSRC theory to the inversion of noisy synthetic and field data using a heuristic algorithm. The multirealization results of the inversion strategy are consistent with the feasibility analysis and demonstrate the potential of the outlined method for practical application.

Uncertainty quantification in the context of seismic imaging is important for interpretinginverted subsurface models and updating reservoir models. The limited illumination, noisydata and poor initial model in the seismic full waveform inversion (FWI) lead to inversionuncertainties. This is particularly true for anisotropic elastic FWI, which suffers from extra parameter trade-off problems. In this work, we address the uncertainty quantificationof anisotropic elastic FWI problem in the framework of Bayesian inference. Specially, weestimate the uncertainties of the subsurface elastic parameters in the Bayesian anisotropicelastic FWI by combining the iterated extended Kalman filter with an explicit representation of the sensitivity matrix with Green’s functions. The sensitivity matrix is based onthe integral equation approach, which is also within the context of nonlinear inverse scattering theory. We give the results of numerical tests with examples for anisotropic elasticmedia. They show that the proposed Bayesian inversion method can provide reasonablereconstructed results for the elastic coefficients of the stiffness tensor and the framework issuitable for accessing the uncertainties.

We have developed linearized complex eikonal equations for orthorhombic media and the analytic solution for homogeneous acoustic media. The linearized formulations are developed from the highly nonlinear complex eikonal equation using a perturbation analysis. The analytic solution is obtained by expanding the complex traveltime in a Taylor series and constructing its coefficients. A critical factor for deriving this analytic expression is to apply the complex point-source method to the analytic expression for the background anisotropic medium. The analytic expressions show that the complex traveltime is a function of the anisotropic parameters, wave velocity, and spatial coordinates. We give the numerical results of the complex traveltime for different anisotropic parameters and initial beam widths. Furthermore, numerical examples are investigated for various coefficients as a function of the coordinates.