Precise modeling of the attenuation parameter Q is important to confirm the robust and diagnostic attenuation characteristics of seismic waveforms in oil and gas reservoirs. For this purpose, the nearly constant Q viscoacoustic wave equation is beneficial. Compared with attenuation models such as standard linear solid (SLS), the wave equation corresponding to the nearly constant Q model contains an explicit Q, which is conducive to inversion and can provide a more effective parameterization. This parameterization facilitates the initial suppression of parameter crosstalk in multiparameter inversion. We derive the adjoint equation containing auxiliary variables, compare the sensitive kernels of different parameterizations, and explain the superiority of the constructed parameterization. The truncated Gauss-Newton (GN-TRN) method is introduced to suppress parameter crosstalk further. The GN-TRN method updates the model parameters by calculating the Hessian vector product and iteratively solving the approximation of the Newton gradient directions. The test results of the theoretical model and field data verify the effectiveness and stability of the inversion method.

Full waveform inversion (FWI) is a powerful means for accurately reconstructing subsurface velocity models at high resolution. Yet it is nevertheless a nonlinear and ill-posed problem. Physics-driven FWI methods employ gradient-based optimization algorithms to minimize the error between the observed seismic data and the synthetically generated seismic data. The solution may converge to a local rather than global minimum. The cycle-skipping problem occurs when the synthetic data exceed a half-wavelength shift relative to the observed data. FWI relies on an accurate initial velocity model to mitigate the cycle-skipping problem. Moreover, due to the increasing size and desired resolution of seismic data, FWI costs a great deal of computational time. To obviate these problems, we present a data-driven FWI scheme based on a deep learning architecture called U-Net. The network consists of the ringed residual unit, which integrates residual propagation and residual feedback. It beneficially achieves correspondence between the seismic data domain and the velocity model domain. The features of the shallow layers are connected with the deep layers by a skip connection to facilitate seismic data spatial information propagation and utilization. They improve inversion accuracy and make the network more generalizable and robust. We utilize the Society of Exploration Geophysicists (SEGs)/European Association of Geoscientists and Engineers (EAGE) overthrust and salt models to verify our proposed method’s impressive performance. The experimental results clearly demonstrate that the proposed method can produce high-quality velocity models. Compared with the conventional physics-informed FWI, it has advantages in both computational time and initial model dependence.

Seismic wavefield forward modeling in anelastic (attenuating) media is a fundamental tool for both data processing and interpretation in modern seismic exploration. We propose a generalized recursive convolution (RC) formula to calculate the temporal convolutions directly, rather than solving many auxiliary partial differential equations of the memory variables when dealing with a viscoelastic medium. The new formula is obtained in terms of the Taylor series expansion and offers approximations of the convolutions to arbitrary order by the number of terms retained. We conduct theoretical and numerical comparisons of the new method with the commonly used memory variable method and other traditional RC methods. The comparisons show that the new method has the highest accuracy of all these RC methods and yields better performance with various stress relaxation times and time steps than the common leapfrog time-stepping scheme to solve the auxiliary partial differential equations of the memory variables. Our numerical examples verify the versatility and feasibility of the new method for viscoacoustic and viscoelastic wave modeling.

The Earth is often characterized by viscoelastic rocks, porous sediments and anisotropic structures. Poro-elasticity with Biot’s theory is considered fundamental to describe the interaction between the deformation of the elastic porous solid and the flow of fluid in the porous structure. The quality factor (Q) in the theory of viscoelasticity relates seismic wave attenuation and dispersion to physical properties of the Earth’s interior, e. g. temperature, stress and composition. However, the constantQwave equation in its time-domain differential form remains difficult to solve when describing the attenuation in an explicitly specifiedQparameter. Here, we introduce the first-and second-order nearly constantQmodels capable of describing the attenuation of the solid skeleton, thereby extending the Biot and Biot-squirt (BISQ) models to poro-viscoelastic media. The bulk and shear moduli of the solid frame are represented by the modified relaxation function. By presenting examples with finite-difference time-domain (FDTD) numerical modeling for seismic wavefields in anisotropic, viscoelastic porous media including transversely isotropic media with a vertical symmetry axis (VTI) and orthorhombic meida, we demonstrate that the extended Biot and BISQ models provide good descriptions of the wave propagation in poro-viscoelastic anisotropic media and can thus help better understand the Earth’s interior.

In this work, a Robust Worst-Case (RWC) estimation is presented for recovering sparse reflectivity series from uniformly quantized seismic signals. First, theOrthogonal Matching Pursuit(OMP) algorithm is applied on a quantized seismic trace. A set of conservative estimates of the reflectivity impulses and a dimension-reduced system are accordingly obtained. Second, the error induced by the quantization is modeled as a systemic uncertainty in a multiplication way. This modeling imposes a greater model uncertainty on the estimated reflectivity impulses with small values and vice versa. Finally, a RWC deconvolution scheme is designed for the dimension-reduced system. Among those roughly estimated impulses from the OMP algorithm, the ones statistically less affected by the quantization process are assigned with higher confidence weights and the ones statistically with higher magnitudes of quantization error are assigned with lower confidence weights. By rescaling the quantization error, the proposed scheme significantly increases the robustness of the solution to the quantization error and hence better improves the visual saliency of seismic signals than that of the OMP algorithm. This scheme is tested on both synthetic and real seismic data and the performance is evaluated by compared to that of the OMP algorithm. The results shows that the new scheme significantly outperforms the OMP algorithm by observing the following two aspects: first, the falsely or overly estimated impulses are significantly suppressed in the RWC estimation compared with that of the OMP algorithm, and second, the RWC estimation exhibits enhanced robustness to the change of the quantization interval.

Radar tomography is a useful technique for mapping the permittivity and conductivity distributions in the shallow subsurface. By exploiting the full radar waveforms, it is possible to improve resolution and, thus, image subwavelength features not resolvable using ray-based approaches. Usually, mere convergence in the data space is the only criterion used to appraise the goodness of a final result, possibly limiting the reliability of the inversion. A better indication of the correctness of an inverted model and its various parts could be obtained by means of a formal model resolution and information content analysis. We present a novel method for computing the sensitivity functions (Jacobian matrix) based on a time-domain adjoint method. Because the new scheme only computes the sensitivity values for the transmitter and receiver combinations that are used, it reduces the number of forward runs with respect to standard brute-force or other virtual-source schemes. The procedure has been implemented by using a standard finite-difference time-domain modeling method. A comparison between cumulative sensitivity (column sum of absolute values of the Jacobian) images, which is sometimes used in geoelectrical studies as a proxy for resolution in practical cases, and formal model resolution images is also presented. We show that the cumulative sensitivity supplies some valuable information about the image, but when possible, formal resolution analyses should be performed. The eigenvalue spectrum of the pseudo-Hessian matrix provides a measure of the information content of an experiment and shows the extent of the unresolved model space.

We have developed a numerical algorithm for computing the magnetic field distribution and the nuclear magnetic resonance (NMR) sensitivity function for arbitrary topography overlying a known 3-D conductivity structure. The magnetic vector potential is split into primary and secondary terms. The primary term is obtained using a thin-wire line integral equation that accounts for arbitrary loop shape and position. It allows the singularity of the source field to be effectively removed. The secondary potential is obtained by solving the second-order partial differential equations on an unstructured tetrahedral mesh using the finite element technique. We validate the results of applying our algorithm against an explicit infinite integral solution for circular loops on a layered earth and against the results of applying a commercial simulation tool. The spatially oscillating NMR sensitivity functions to hydrogen protons (i.e., unbounded water molecules) in the sub-surface are computed on a refined unstructured grid. We apply the numerical algorithm to a number of synthetic examples in surface NMR tomography of hydrological relevance.

We have developed a new full-waveform groundpenetrating radar (GPR) multicomponent inversion scheme for imaging the shallow subsurface using arbitrary recording configurations. It yields significantly higher resolution images than conventional tomographic techniques based on first-arrival times and pulse amplitudes. The inversion is formulated as a nonlinear least squares problem in which the misfit between observed and modeled data is minimized. The full-waveform modeling is implemented by means of a finite-difference time-domain solution of Maxwell's equations. We derive here an iterative gradient method in which the steepest descent direction, used to update iteratively the permittivity and conductivity distributions in an optimal way, is found by cross-correlating the forward vector wavefield and the backward-propagated vectorial residual wavefield. The formulation of the solution is given in a very general, albeit compact and elegant, fashion. Each iteration step of our inversion scheme requires several calculations of propagating wavefields. Novel features of the scheme compared to previous full-waveform GPR inversions are as follows: 1) The permittivity and conductivity distributions are updated simultaneously (rather than consecutively) at each iterative step using improved gradient and step length formulations; 2) the scheme is able to exploit the full vector wavefield; and 3) various data sets/survey types (e.g., crosshole and borehole-to-surface) can be individually or jointly inverted. Several synthetic examples involving both homogeneous and layered stochastic background models with embedded anomalous inclusions demonstrate the superiority of the new scheme over previous approaches.

This article addresses the comments of Sakamato et al. on Greenhalgh et al.'s "Modeling and migration of 2-D georadar data: A stationary phase approach". Here, Greenhalgh et al. emphasizes the comprehensiveness of their approach over Sakamato et al.'s SEABED algorithm.

This paper presents the basic kinematic and dynamic imaging and migration equations for zero-offset two-dimensional georadar profiling. The kinematic equations are derived from simple considerations of spatial impulse responses and a generating function. The dynamic equations follow from a multidimensional stationary phase approximation to the infinite spectral integrals. They show how the radar signal (amplitude and phase) depends on the shape and curvature of the reflector. The imaging equations are evaluated for the special cases of a point scatterer, a continuous reflector, and a terminating reflector. A general formula is developed by which to migrate an arbitrary shaped event of variable amplitude on the georadar section