Ru-Shan Wu is a special scientist who has made tremendous contributions in many fields of geophysics that he has worked on and has had strong impact on several generations of scientists throughout the world. He is appreciated for his contributions to science and for his generous, outgoing and happy personality that make him a role model for those who follow in his footsteps. Ru-Shan has sustained high-quality scientific output over the course of his long career. His work is widely recognized for its value and creativity. He has keen insight into the underlying physics of wave propagation and he continues to successfully transform that insight into fundamental advances in our understanding of wave propagation and imaging.

Ru-Shan Wu has made seminal contributions in many research areas in geophysics, such as seismic-wave propagation, scattering, imaging, and inversion. We highlight some of his research in holography imaging, diffraction tomography, seismic-wave scattering and its applications to studying Earth's heterogeneity, one-way wave propagation and one-return wave modeling, beamlet and dreamlet applications, strong non-linear full-waveform inversion, and direct envelop inversion.

Strong-scattering inversion or the inverse problem for strong scattering has different physical-mathematical foundations from the weak-scattering case. Seismic inversion based on wave equation for strong scattering cannot be directly solved by Newton's local optimization method which is based on weak-nonlinear assumption. Here I try to illustrate the connection between the Schrödinger inverse scattering (inverse problem for Schrödinger equation) by GLM (Gel'fand-Levitan-Marchenko) theory and the direct envelope inversion (DEI) using reflection data. The difference between wave equation and Schrödinger equation is that the latter has a potential independent of frequency while the former has a frequency-square dependency in the potential. I also point out that the traditional GLM equation for potential inversion can only recover the high-wavenumber components of impedance profile. I propose to use the Schrödinger impedance equation for direct impedance inversion and introduce a singular impedance function which also corresponds to a singular potential for the reconstruction of impedance profile, including discontinuities and long-wavelength velocity structure. I will review the GLM theory and its application to impedance inversion including some numerical examples. Then I analyze the recently developed multiscale direct envelope inversion (MS-DEI) and its connection to the inverse Schrödinger scattering. It is conceivable that the combination of strong-scattering inversion (inverse Schrödinger scattering) and weak-scattering inversion (local optimization based inversion) may create some inversion methods working for a whole range of inversion problems in geophysical exploration.

Fluctuations in log-amplitude and travel time of teleseismic P waves recorded by the EarthScope USArray are used to invert for the heterogeneity spectrum of P-wave velocity in a 1000 km thick region of the upper mantle beneath the array. These fluctuations are used to form coherence functions. Best fits to joint transverse coherence functions require a depth dependent heterogeneity spectrum, with peaks in narrow depth ranges. These peaks agree well with peaks predicted for the temperature derivative of seismic velocity from models of the chemistry and phase of silicate mineral assemblages appropriate for the upper mantle, correlating with the depths of phase changes. The results show that chemistry and phase act in concert with lateral and depth variations in temperature to produce the observed heterogeneity in seismic velocities in the upper mantle at spatial scales from 50 to 300 km.

In the field of geophysics, although the first-order Rytov approximation is widely used, the higher-order approximation is seldom discussed. From both theoretical analysis and numerical tests, the accumulated phase error introduced in the first-order Rytov approximation cannot be neglected in the presence of strong velocity perturbation. In this paper, we are focused on improving the phase accuracy of forward scattered wavefield, especially for the large-scale and strong velocity perturbation case. We develop an equivalent source method which can update the imaginary part of the complex phase iteratively, and the higher-order scattered wavefield can be approximated by multiplying the incident wavefield by the exponent of the imaginary part of the complex phase. Although the convergence of the proposed method has not been proved mathematically, numerical examples demonstrate that our method can produce an improved accuracy for traveltime (phase) prediction, even for strong perturbation media. However, due to the neglect of the real part of the complex phase, the amplitude change of the scattered wavefield cannot be recovered. Furthermore, in the presence of multi-arrivals phenomenon, the equivalent scattering source should be handled carefully due to the multi-directions of the wavefield. Further investigations should be done to improve the applicability of the proposed method.

Seismic events have limited time duration, vary with space/traveltime and interact with the local subsurface medium during propagation. Partitioning is a valuable strategy for nonstationary seismic data analysis, processing and wave propagation. It has the potential for sparse data representation, flexible data operation and highly accurate local wave propagation. Various local transforms are powerful tools for seismic data segmentation and representation. In this paper, a detailed description of a multi-dimensional local harmonic transformed domain wave propagation and imaging method is given. Using a tensor product of a Local Exponential Frame (LEF) vector as the time-frequency atom (a drumbeat) and a Local Cosine Basis (LCB) function as the space-wavenumber atom (a beamlet), we construct a time-frequency-space-wavenumber local atom-dreamlet, which is a combination of drumbeat and beamlet. The dreamlet atoms have limited spatial extension and temporal duration and constitute a complete set of frames, termed as dreamlet frames, to decompose and represent the wavefield. The dreamlet transform first partitions the wavefields using time-space supporting functions and then the data in each time-space blocks is represented by local harmonic bases. The transformed wavefield is downward-continued by the dreamlet propagator, which is the dreamlet atom evolution weightings deduced from the phase-shift one-way propagator. The dreamlet imaging method is formulated with a local background propagator for large-scale medium propagation and combined with a local phase-screen correction for small-scale perturbations. The features of dreamlet migration and imaging include sparse seismic data representation, accurate wave propagation and the flexibility of localized time operations during migration. Numerical tests using Sigsbee 2A synthetic data set and real marine seismic data demonstrate the validity and accuracy of this method. With time-domain localization being involved, the dreamlet method can also be applied effectively to target-oriented migration and imaging.

Converting subsurface offset domain common image gathers (ODCIGs) to angle domain common image gathers (ADCIGs) through a Radon Transform (RT) in either the spatial or wavenumber domain is efficient and valid except for the distortion of both frequency spectrum and amplitude versus angle (AVA) effect. This paper presents two modifications to the existing method to keep the frequency spectrum of the resultant ADCIGs the same as the input data and to preserve the relative amplitudes. The spectrum invariance is achieved by replacing the conventional RT or slant slack by an invertible RT. Amplitude preservation is obtained by applying an amplitude correction factor in the angle domain. Tests on both synthetic and field datasets validate the accuracy of these modifications.

The stability and efficiency, especially the stability, are generally concerned issues in Q compensated reverse time migration (Q-RTM). The instability occurs because of the exponentially boosted high frequency ambient noise during the forward or backward seismic wavefield propagation. The regularization and low-pass filtering methods are two effective strategies to control the instability of the wave propagation in Q-RTM. However, the regularization parameters are determined experimentally, and the wavefield cannot be recovered accurately. The low-pass filtering method cannot balance the selection of cutoff frequency for varying Q values, and may damage the effective signals, especially when the signal-to-noise ratio (SNR) of the seismic data is low, the Q-RTM will be a highly unstable process. In order to achieve the purpose of stability, the selection of cutoff frequency will be small enough, which can cause great damage to the effective high frequency signals. In this paper, we present a stable Q-RTM algorithm based on the excitation amplitude imaging condition, which can compensate both the amplitude attenuation and phase dispersion. Unlike the existing Q-RTM algorithms enlarging the amplitude, the exponentially attenuated seismic wavefield will be used during both the forward and backward wavefield propagation of Q-RTM. Therefore, the new Q-RTM algorithm is relative stable, even for the low SNR seismic data. In order to show the accuracy and stability of our stable Q-RTM algorithm clearly, an example based on Graben model will be illustrated. Then, a realistic BP gas chimney model further demonstrates that the proposed method enjoys good stability and anti-noise performance compared with the traditional Q-RTM with amplitude amplification. Compare the Q-RTM images of these two models to the reference images obtained by the acoustic RTM with acoustic seismic data, the new Q-RTM results match the reference images quite well. The proposed method is also tested using a field seismic data, the result shows the effectiveness of our proposed method.

In land seismic exploration, strong near-surface heterogeneities can cause serious problems in seismic data acquisition and the quality of depth imaging. By introducing random velocity models to simulate velocity fluctuations in the near-surface layer and using the point spread function to characterize image quality, we examine how the scattering generated in near-surface heterogeneities can affect the subsurface image. In addition to the commonly known scattering noises which lower the signal to noise ratio in seismic data, our results also reveal that intermediate scale heterogeneities generate forward scattering which forms phase or travel time fluctuations. Due to intermediate-scale uncertainty in the shallow part of the migration velocity model, these phase changes are carried to the target by the extrapolated wavefields, breaking the zero phase image condition at the image point. This is a primary reason for deteriorated image quality in regions with strong near-surface scattering. If this intermediate-scale information can be obtained and built into the migration velocity, the subsurface image quality can be largely improved. These results can be the basis for further numerical investigations and field experiments. The proposed analysis method can also be used to evaluate other potential methods for dealing with near-surface scattering.

Staining algorithms based on two-way wave equation migration methods have been applied to improve the signal-to-noise ratio (SNR) of poorly illuminated structures such as those in subsalt zones. In regular staining algorithms, when a source wavefield reaches the stained area that is associated with the target structures, a new wavefield called stained wavefield is excited, and this stained wavefield forward extrapolates synchronously with the real source wavefield. The forward-extrapolated stained and real source wavefields are cross-correlated with the backward-extrapolated receiver wavefield, and we obtain the stained and the real reverse time migration (RTM) images. The staining algorithms for RTM can suppress the noise of non-target regions and obtain high SNR images of the target structures. Whereas RTM methods are limited by the low computational efficiency and SNR, by contrast, one-way wave equation migration (OWEM) methods have the advantages of high efficiency and no interference from multiples. Thus, we developed a new staining method based on the generalised screen propagator (GSP) as a case of OWEM methods for subsalt imaging. Furthermore, a new stained wavefield called stained receiver wavefield is proposed here, forming two new staining strategies for seismic imaging, in which forward-propagated source and backward-propagated receiver wavefields can be conveniently selected to be stained at the stained area. Numerical experiments demonstrated that this staining GSP method is more effective in improving the SNR of subsalt structures compared to conventional GSP migration and RTM methods; moreover, these new staining strategies as applied to the OWEM methods can greatly improve the SNR of weakly illuminated structures in subsalt zones, in comparison with regular staining algorithms for one-way methods.

In the contrast source inversion (CSI) method, the contrast sources (equivalent scattering sources) and the contrast (parameter perturbation) are iteratively reconstructed by an alternating optimization scheme. Traditionally integral equation CSI method is formulated for transmission tomography using analytic Green's function in homogeneous background. To extend the method to the case of reflection seismology, in this paper, we use WKBJ method to compute the Green's function of depth dependent background media and the solving method of equations to initialize the contrast source of different frequencies, resulting in an efficient method to invert multi-frequency reflection seismic data – multi-frequency contrast source inversion method (MFCSI). Numerical results for the Marmousi model show that MFCSI method can obtain good results even when low frequency data are missing, in which case the conventional FWI fails.

In full waveform inversion (FWI) high-resolution subsurface model parameters are sought. FWI is normally treated as a nonlinear least-squares inverse problem, in which the minimum of the corresponding misfit function is found by updating the model parameters. When multiple elastic or acoustic properties are solved for, simple gradient methods tend to confuse parameter classes. This is referred to as parameter cross-talk; it leads to incorrect model solutions, poor convergence and strong dependence on the scaling of the different parameter types. Determining step lengths in a subspace domain, rather than directly in terms of gradients of different parameters, is a potentially valuable approach to address this problem. The particular subspace used can be defined over a span of different sets of data or different parameter classes, provided it involves a small number of vectors compared to those contained in the whole model space. In a subspace method, the basis vectors are defined first, and a local minimum is found in the space spanned by these. We examine the application of the subspace method within acoustic FWI in determining simultaneously updates for velocity and density. We first discuss the choice of basis vectors to construct the spanned space, from linear updates by distinguishing only the contributions of different parameter classes towards nonlinear updates by adding the contributions of higher-order perturbations of each parameter class. The numerical character of FWI solutions generated via subspace methods involving different basis vectors is then analyzed and compared with traditional FWI methods. The subspace methods can provide better reconstructions of the model, especially for the velocity, as well as improved convergence rates, while the computational costs are still comparable with the traditional FWI methods.

Full waveform inversion of time-lapse seismic data can be used as a means of estimating the reservoir changes due to the production. Since the repeated computations for the monitor surveys lead to a large computational cost, time-lapse full waveform inversion is still considered to be a challenging task. To address this problem, we present an efficient target-oriented inversion scheme for time-lapse seismic data using an integral equation formulation with Gaussian beam based Green's function approach. The proposed time-lapse approach allows one to perform a local inversion within a small region of interest (e.g. a reservoir under production) for the monitor survey. We have verified that the T-matrix approach is indeed naturally target-oriented, which was mentioned by Jakobsen and Ursin [24] and allows one to reduce the computational cost of time-lapse inversion by focusing the inversion on the target-area only. This method is based on a new version of the distorted Born iterative T-matrix inverse scattering method. The Gaussian beam and T-matrix are used in this approach to perform the wavefield computation for the time-lapse inversion in the baseline model from the survey surface to the target region. We have provided target-oriented inversion results of the synthetic time-lapse waveform data, which shows that the proposed scheme reduces the computational cost significantly.

Recorded seismic data are sensitive to the Earth's elastic properties, and thus, they carry information of such properties in their waveforms. The sensitivity of such waveforms to the properties is nonlinear causing all kinds of difficulties to the inversion of such properties. Inverting directly for the components forming the wave equation, which includes the wave equation operator (or its perturbation), and the wavefield, as independent parameters enhances the convexity of the inverse problem. The optimization in this case is provided by an objective function that maximizes the data fitting and the wave equation fidelity, simultaneously. To enhance the practicality and efficiency of the optimization, I recast the velocity perturbations as secondary sources in a modified source function, and invert for the wavefield and the modified source function, as independent parameters. The optimization in this case corresponds to a linear problem. The inverted functions can be used directly to extract the velocity perturbation. Unlike gradient methods, this optimization problem is free of the Born approximation limitations in the update, including single scattering and cross talk that may arise for example in the case of multi sources. These specific features are shown for a simple synthetic example, as well as the Marmousi model.

We generalize the existing distorted Born iterative T-matrix (DBIT) method to seismic full-waveform inversion (FWI) based on the scalar wave equation, so that it can be used for seismic FWI in arbitrary anisotropic elastic media with variable mass densities and elastic stiffness tensors. The elastodynamic wave equation for an arbitrary anisotropic heterogeneous medium is represented by an integral equation of the Lippmann-Schwinger type, with a 9-dimensional wave state (displacement-strain) vector. We solve this higher-dimensional Lippmann-Schwinger equation using a transition operator formalism used in quantum scattering theory. This allows for domain decomposition and novel variational estimates. The tensorial nonlinear inverse scattering problem is solved iteratively by using an expression for the Fréchet derivatives of the scattered wavefield with respect to elastic stiffness tensor fields in terms of modified Green's functions and wave state vectors that are updated after each iteration. Since the generalized DBIT method is consistent with the Gauss-Newton method, it incorporates approximate Hessian information that is essential for the reduction of multi-parameter cross-talk effects. The DBIT method is implemented efficiently using a variant of the Levenberg-Marquard method, with adaptive selection of the regularization parameter after each iteration. In a series of numerical experiments based on synthetic waveform data for transversely isotropic media with vertical symmetry axes, we obtained a very good match between the true and inverted models when using the traditional Voigt parameterization. This suggests that the effects of cross-talk can be sufficiently reduced by the incorporation of Hessian information and the use of suitable regularization methods. Since the generalized DBIT method for FWI in anisotropic elastic media is naturally target-oriented, it may be particularly suitable for applications to seismic reservoir characterization and monitoring. However, the theory and method presented here is general.

Full-waveform inversion is a promising tool to produce accurate and high-resolution subsurface models. Conventional full-waveform inversion requires an accurate estimation of the source wavelet, and its computational cost is high. We develop a novel source-independent full-waveform inversion method using a hybrid time- and frequency-domain scheme to avoid the requirement of source wavelet estimation and to reduce the computational cost. We employ an amplitude-semblance objective function to not only effectively remove the source wavelet effect on full-waveform inversion, but also to eliminate the impact of the inconsistency of source wavelets among different shots gathers on full-waveform inversion. To reduce the high computational cost of full-waveform inversion in the time domain, we implement our new algorithm using a hybrid time- and frequency-domain approach. The forward and backward wave propagation operations are conducted in the time domain, while the frequency-domain wavefields are obtained during modeling using the discrete-time Fourier transform. The inversion process is conducted in the frequency domain for selected frequencies. We verify our method using synthetic seismic data for the Marmousi model. The results demonstrate that our novel source-independent full-waveform inversion produces accurate velocity models even if the source signature is incorrect. In addition, our method can significantly reduce the computational time using the hybrid time- and frequency-domain approach compared to the conventional full-waveform inversion in the time domain.

The Direct Waveform Inversion (DWI) is a recently proposed waveform inversion idea that has the potential to simultaneously address several existing challenges in many full waveform inversion (FWI) schemes. A key ingredient in DWI is the explicit use of the time-space causality property of the wavefield in the inversion which allows us to convert the global nonlinear optimization problem in FWI, without information loss, into local linear inversions that can be readily solved. DWI is a recursive scheme which sequentially inverts for the subsurface model in a shallow-to-deep fashion. Therefore, there is no need for a global initial velocity model to implement DWI. DWI is unconditionally convergent when the reflection traveltime from the boundary of inverted model is beyond the finite recording time in seismic data. In order for DWI to work, DWI must use the full seismic wavefield including interbed and free surface multiples and it combines seismic migration and velocity model inversion into one process. We illustrate the concepts in DWI using 1D and 2D models.

The goals of this study were to examine factors influencing Q inversion and to provide references for practical application. Three different methods for inverting Q values with VSP data were explored, including centroid frequency shift (CFS), spectral ratio (SR), and amplitude attenuation (AA). Comparison between the CFS and the other two methods was conducted on frequency band widths and low attenuation, wavefield components, interface interference, and thin layers. Results from several sets of VSP modeling data indicated that the CFS method is more stable and accurate for dealing with thin and high Q layers. Frequency band width, especially the presence of high frequencies, influences the inversion effect of all three methods. The wider the band, the better the results. Q inversion from downgoing wavefield was very similar to that of the upgoing wavefield. The CFS method had fewer outliers or skip values from the full wavefield than the other two methods. Moreover, the applications to Q inversion for the set of field VSP data demonstrated that the Q curves from the CFS method coincided with the geological interpretations better than the Q curves of the other methods. Meanwhile, inverse Q filtering shifted the frequency component from 25 Hz to 35 Hz. The results demonstrated that the Q curve is more sensitive to geological horizons than velocity.

Numerically solving 3D seismic wave equations is a key requirement for forward modeling and inversion. Here, we propose a weighted Runge-Kutta discontinuous Galerkin (WRKDG) method for 3D acoustic and elastic wave-field modeling. For this method, the second-order seismic wave equations in 3D heterogeneous anisotropic media are transformed into a first-order hyperbolic system, and then we use a discontinuous Galerkin (DG) solver based on numerical-flux formulations for spatial discretization. The time discretization is based on an implicit diagonal Runge-Kutta (RK) method and an explicit iterative technique, which avoids solving a large-scale system of linear equations. In the iterative process, we introduce a weighting factor. We investigate the numerical stability criteria of the 3D method in detail for linear and quadratic spatial basis functions. We also present a 3D analysis of numerical dispersion for the full discrete approximation of acoustic equation, which demonstrates that the WRKDG method can efficiently suppress numerical dispersion on coarse grids. Numerical results for several different 3D models including homogeneous and heterogeneous media with isotropic and anisotropic cases show that the 3D WRKDG method can effectively suppress numerical dispersion and provide accurate wave-field information on coarse mesh.

Numerical solution of time-lapse seismic monitoring problems can be challenging due to the presence of finely layered reservoirs. Repetitive wave modeling using fine layered meshes also adds more computational cost. Conventional approaches such as finite difference and finite element methods may be prohibitively expensive if the whole domain is discretized with the cells corresponding to the grid in the reservoir subdomain. A common approach in this case is to use homogenization techniques to upscale properties of subsurface media and assign the background properties to coarser grid; however, inappropriate application of upscaling might result in a distortion of the model, which hinders accurate monitoring of the fluid change in subsurface. In this work, we instead investigate capabilities of a multiscale method that can deal with fine scale heterogeneities of the reservoir layer and more coarsely meshed rock properties in the surrounding domains in the same fashion. To address the 3-D wave problems, we also demonstrate how the multiscale wave modeling technique can detect the changes caused by fluid movement while the hydrocarbon production activity proceeds.

Evaluation of possible leakage pathways of CO₂ injected into geological formations for storage is essential for successful Carbon Capture and Storage (CCS). A channel in the borehole cement, which secures the borehole casing to the formation, may allow CO₂ to escape. Risk assessment and remediation decisions about the presence of such channels depend on channel parameters: radial position $r$ from the center of the borehole; channel thickness $d$; azimuthal position $φ$ of the channel; and azimuthal extent $θ$ of the channel. Current state-of-the-art cement-bond logging technology, which uses only the first arrival at a centralized borehole receiver, can diagnose limited details about CO₂ leak channels. To accurately characterize the possible leak paths in the cement, we use a 3-dimensional finite-difference method to investigate the use of the abundant data collected by a modernized monopole sonic tool that contains an array of azimuthally distributed receivers. We also investigate how to improve the tool design to acquire even more useful information. For cases where borehole fluid is either water or supercritical CO₂, we investigate various receiver geometries, multimodal analyses of multi-frequency data to discover the type of logging tool that provides the best information for CCS management. We find that an appropriate choice of wave modes, source frequencies, source polarities, and receiver locations and offsets provides sensitivity to $d$, $φ$, $θ$. The amplitude of the first arrival from a monopole source is sensitive to $θ$. Amplitudes at receivers at different azimuths are sensitive to $φ$. The slow Stoneley mode (ST2) velocity is sensitive to $d$, but ST2 is not easy to pick when $θ$ and $d$ are small. Further improvement is necessary to provide comprehensive information about possible flow channels in casing cement.

Reliable subsurface time-lapse seismic monitoring is crucial for many geophysical applications, such as enhanced geothermal system characterization, geologic carbon utilization and storage, and conventional and unconventional oil/gas reservoir characterization, etc. We develop an elastic-wave sensitivity propagation method for optimal design of cost-effective time-lapse seismic surveys considering the fact that most of subsurface geologic layers and fractured reservoirs are anisotropic instead of isotropic. For anisotropic media, we define monitoring criteria using qP- and qS-wave sensitivity energies after decomposing qP- and qS-wave components from the total elastic-wave sensitivity wavefield using a hybrid time- and frequency-domain approach. Geophones should therefore be placed at locations with significant qP- and qS-wave sensitivity energies for cost-effective time-lapse seismic monitoring in an anisotropic geology setting. Our numerical modeling results for a modified anisotropic Hess model demonstrate that, compared with the isotropic case, subsurface anisotropy changes the spatial distributions of elastic-wave sensitivity energies. Consequently, it is necessary to consider subsurface anisotropies when designing the spatial distribution of geophones for cost-effective time-lapse seismic monitoring. This finding suggests that it is essential to use our new anisotropic elastic-wave sensitivity modeling method for optimal design of time-lapse seismic surveys to reliably monitor the changes in subsurface reservoirs, fracture zones or target monitoring regions.

Removing internal multiples remains an important but challenging problem in seismic processing. The generalized Estimation of Primaries by Sparsity Inversion (EPSI) method minimizes data residuals between the calculated and observed waveform using the sparse constraint of primary impulse responses to predict multiples and remove them directly, instead of using the conventional adaptive subtraction method. Even though the generalized EPSI method provides a good estimate of the primaries and multiples when they overlap, it is limited by intensive computational cost.
In this paper, we introduce two strategies to improve computational efficiency. First, the interface-controlled strategy is introduced by only selecting high-amplitude primary responses related to the interfaces with strong impedance contrasts to estimate multiples. The computational time is approximately proportional to the number of involved reflectors and usually, most of the internal multiple energy in the data is only related to a few strong reflectors. Therefore the modified method can remove most of the internal multiples in fewer computations than in the generalized EPSI, which loops through all the interfaces. Next, an approximate formula for estimating primary impulse responses is proposed by neglecting a computationally intensive term which corresponds to the primary responses estimated from internal multiples. According to our analyses and experiments, in most cases, the contribution of this term is negligible because the internal multiples are weak. Therefore, the computational efficiency can be improved without significantly losing quality when estimating most primaries and multiples.
In order to demonstrate this, multiple elimination of a two-layered simple data and the Pluto data are implemented. We find that the modified method can yield reliable results that require fewer computations. The improvements of the modified method may encourage the use of the generalized EPSI method in industry.

Recently, attenuation of surface-related multiples is implemented by a large-scale sparsity-promoting inversion where the primaries are iteratively estimated without a subtraction process, which is called estimation of primaries by sparse inversion (EPSI). By inverting for surface-free impulse responses, EPSI simultaneously updates the primaries and multiples, both of which contribute to explaining the input data, and therefore promote the global convergence gradually. However, one of the major concerns of EPSI may lie in its high computational cost. In this paper, based on the same gradient-descent framework with EPSI, we develop a computationally cost-effective primary estimation approach in which a newly defined parameterization of primary-multiple model is adopted and an efficiently defined analytical step-length is developed. The developed approach can yield a better primary estimation at less computational cost as compared to EPSI, which is verified by two synthetic datasets in numerical examples. Moreover, we apply this approach to a shallow-water field dataset and achieve a desirable performance.

This study is aimed at improving multiple adaptive subtraction. We propose a modified pseudomulti-channel matching method based on the Huber norm, to adjust the matching differences on frequency and phase between the predicted multiples and original data. The second-order derivative of the predicted multiples is utilized to replace the derivative of its Hilbert transform. Due to the additional frequency term, this method can enhance the high-frequency component. We introduce 180〫 phase rotation of the multiple channels, which can decrease phase differences. The Huber norm interpolates between smooth L2 norm treatment of small residuals and robust L1 norm treatment of large residuals. This method can eliminate the restriction of large value conditions from the L2 norm and weaken the condition of orthogonality from the L1 norm. The applications of the Pluto and Delft models shows that compared with pseudomulti-channel matching filter, the main frequency is increased from 36 Hz to 38 Hz, and the primary reflection wave is more concentrated. The practical application of field data verifies the effectiveness of the proposed method.

Effective medium methods for the attribution of micro-structures to macro elastic properties of shales are important for the prediction of sweet spots in the shale-gas production. With X-ray micro-computed tomography (XMCT), the micro-structures of shale core samples from Longmaxi Formation are visualized and characterized by 3D digital images. As an efficient alternative to conventional effective medium methods for estimating elastic properties, we propose a consistent workflow of lattice spring modeling (LSM) to emulate the digital cores using three types of lattices. Particular attention is paid to investigate the effective Young's moduli, Poisson's ratios, and preferred orientations, by uniaxial compression tests along two directions. Within elastic deformation, the impact of lattice arrangements on the anisotropy is even more than those of stress disturbances and micro-structural features. Compared with analytical approximations and theoretical predictions, the LSM numerical scheme shows general applicability for heterogeneous porous rocks.