Wave Equation Imaging Research Articles

Assessing the reliability of local earthquake tomography for crustal imaging: 30 yr of records in the Western Alps as a case study

SUMMARY Local earthquake tomography (LET) is a popular method for inverting arrival time picks of local-regional earthquakes for P- and S-wave velocities and hypocentre parameters in seismically active regions. This popularity is due to some robust and well-documented open-source codes that are sometimes used as black boxes. The availability of a very complete time-pick database on the Western Alps gives us a chance to thoroughly investigate the influence of the numerous processes and parameters involved when applying LET to the Western Alps or similar targets. From a subset of high-quality manual picks (1989–2014), we compute preliminary P and S velocity models that are used to predict arrival times for later events and allow the selected fusion of picks downloaded from different seismological agencies for a consistent 33-yr database (1989–2021). Different model reconstructions are proposed by selecting different subsets of the arrival time data set. Aside data input into the LET, influence of initial stratified or 3-D velocity models and related initial earthquake locations is investigated together with grid discretization, Laplacian smoothing and damping parameters in the standard penalty approach commonly used by LET codes. These numerical parameters account for the expected limited resolution of seismic waves due to their finite-frequency content. Parameter selection is handled by the user, whereas frequency impact is only implicit in the onset picks. The earthquake distribution allows a reconstruction down to 40-km depth over an area of ∼500 × 500 km2. Robust features such as the high-velocity Ivrea body anomaly, and a deep low-velocity anomaly associated with crustal thickening underneath the mountain belt survive whatever the tomography strategy and parameters. Finally, a comparison with previous LET reconstructions suggests that finite-frequency content be explicitly incorporated through wave equation tomography to improve spatial resolution. This would fully exploit observables collected from seismograms, albeit with a significant increase in computer costs.

An Envelope Travel-Time Objective Function for Reducing Source–Velocity Trade-Offs in Wave-Equation Tomography

In conventional cross-correlation (CC)-based wave-equation travel-time tomography, wrong source wavelets can result in inaccurate velocity inversion results, which is known as the source–velocity trade-off. In this study, an envelope travel-time objective function is developed for wave-equation tomography to alleviate the non-uniqueness and uncertainty due to wrong source wavelets. The envelope of a seismic signal helps reduce the waveform fluctuations/distortions caused by variations of the source time function. We show that for two seismic signals generated with different source wavelets, the travel-time shift calculated by cross-correlation of their envelopes is more accurate compared to that obtained by directly cross-correlating their waveforms. Then, the CC-based envelope travel-time (ET) objective function is introduced for wave-equation tomography. A new adjoint source has also been derived to calculate the sensitivity kernels. In the numerical inversion experiments, a synthetic example with cross-well survey is first given to show that compared to the traditional CC travel-time objective function, the ET objective function is relatively insensitive to source wavelet variations and can reconstruct the elastic velocity structures more reliably. Finally, the effectiveness and advantages of our method are verified by inversion of early arrivals in a practical seismic survey for recovering near-surface velocity structures.

Quantitative Error Analysis for the Least-Squares Imaging

As oil and gas exploration moves towards complicated geological environments, high-resolution and true-amplitude seismic imaging becomes increasingly important for detecting and evaluating hydrocarbon reservoirs. Traditional ray-based and wave-equation imaging methods can be considered as the adjoint operator of seismic forward modeling, which are difficult to produce high-quality images in complicated structures because of limited frequency band, unbalanced illumination and irregular acquisition. Least-squares migration (LSM) generates an inverse solution for subsurface reflectivity model with high image resolution and balanced amplitudes. Previous studies on LSM mainly focused on the developments of theoretical and practical strategies, but few on error and uncertainty analysis. We present a quantitative analysis method to evaluate the errors of LSM results. The &#x03D5;_<i>data</i> and &#x03C8;_<i>data</i> functions are first computed based on the local similarity and <i>L</i>₂-norm misfit between observed and synthetic data. They are used as data-domain kinematic and dynamic errors, respectively. Then, these local functions are mapped to the subsurface using a Kirchhoff-integral relation to calculate image-domain errors. Numerical examples for synthetic and field data demonstrate that as the iteration of LSM increases, the total kinematic and dynamic errors are reduced, and they vary in different locations. For low signal-to-noise-ratio field data, LSM might enlarge image errors at large iterations because of the overfitting issue, and proper regularization is very important to facilitate the convergence of LSM to a good solution.

Reflection angle/azimuth-dependent least-squares reverse time migration

Seismic images under complex overburdens such as salt are strongly affected by illumination variations due to overburden velocity variations and imperfect acquisition geometries, making it difficult to obtain reliable image amplitudes. Least-squares reverse time migration (LSRTM) addresses these issues by formulating full wave-equation imaging as a linear inverse problem and solving for a reflectivity model that explains the recorded seismic data. Because subsurface reflection coefficients depend on the incident angle, and possibly on the azimuth, quantitative interpretation under complex overburdens requires LSRTM with output in terms of image gathers, e.g., as a function of the reflection angle or angle and azimuth. We have developed a reflection angle- or angle/azimuth-dependent LSRTM method aimed at obtaining physically meaningful image amplitudes interpretable in terms of angle- or angle/azimuth-dependent reflection coefficients. The method is formulated as a linear inverse problem solved iteratively with the conjugate gradient method. It requires an adjoint pair of linear operators for reflection angle/azimuth-dependent migration and demigration based on full wave-equation propagation. We implement these operators in an efficient way by using a mapping approach between migrated shot gathers and subsurface reflection angle/azimuth gathers. To accelerate convergence of the iterative inversion, we apply image-domain preconditioning operators computed from a single de-remigration step. An angle continuity constraint and a structural dip constraint, implemented via shaping regularization, are used to stabilize the solution in the presence of limited illumination and to control the effects of coherent noise. We examine the method on a synthetic data example and on a wide-azimuth streamer data set from the Gulf of Mexico, where we find that angle/azimuth-dependent LSRTM can achieve significant uplift in subsalt image quality, with overburden- and acquisition-related illumination variation effects on angle/azimuth-dependent image amplitudes largely removed.

Estimation of micro-earthquake source locations based on full adjoint P and S wavefield imaging

SUMMARY Locating micro-earthquakes with high resolution and accuracy is a challenge for traveltime inversion, which has uncertainty on the order of a Fresnel zone (many wavelengths). We develop a wave-equation imaging method to increase resolution and reduce location errors to less than a wavelength, but requires very densely deployed receiver arrays with wide aperture and considerable computational cost. Instead of using acoustic data or direct P wave arrivals only, we use elastic multicomponent data and present a new method that uses the full P and S adjoint wavefields to image the microseismic source locations. We separate the P and S waves from the data, and extrapolate the P and S wavefields of each receiver subarray by solving the P and S adjoint wave equations in parallel. We formulate three source imaging conditions by multiplying over subarrays the adjoint P wavefield (IP), S wavefield (IS) and cross-correlated P and S wavefields (IPS). We perform numerical experiments on the highly realistic SEG SEAM4D reservoir model using surface acquisition array geometries. Results for 2-D and 3-D microseismic source estimations show clean images without noisy artefacts at shallow depths. In particular, IPS provides the highest resolution source location image, while IP is limited by the P wavelength and IS is influenced by small coda artefacts. The major-axis alignment and resolution of the source location image are determined by the hypocentral location with respect to the receiver array and illumination-angle coverage, respectively. We discuss the impacts of S-wave attenuation and frequency bandwidth on the source location images. Noise tests indicate that the imaging results are relatively insensitive to ambient noise, as is observed for the surface monitoring data. Using smoothed velocity models, the imaging results are similar to the results using the true realistically heterogeneous velocity model. The 90 per cent confidence ellipse of the source location due to Gaussian-distributed velocity errors shows a larger depth error as the source becomes deeper, while the horizontal error does not change as much.

Low-frequency, long-offset elastic waveform inversion in the context of velocity model building

Classic imaging approaches consist of splitting the earth into background and reflectivity models. When justified, this separation of scale is quite powerful, although this approach relies on some smoothness and weak contrast assumptions. This approach allows for the imaging methods to be based on acoustic wave propagation after having identified the compressional waves through picking or signal processing. Over the past years, wave-equation tomography and waveform inversion approaches have become routine, complementing the classic approaches to derive background models. They do not rely on high-frequency picks, unlike ray-based traveltime tomography, but on low-frequency cross-correlation to define time shifts and on waveform matching. In the presence of large earth parameter contrasts, time shifts and waveforms of compressional waves may depend on elastic parameters when interferences occur within the Fresnel zones. This challenges the recovery of the background model under an acoustic assumption with low-frequency data. Accounting for an elastic propagation in waveform inversion, even in the context of model building, could help to reduce the artifacts seen in acoustic results. A synthetic and a real data example are presented to illustrate the potential benefit of using an elastic waveform inversion approach when inverting long-offset, low-frequency seismic data.

PS energy imaging condition for microseismic data — Part 1: Theory and applications in 3D isotropic media

Accurately estimating event locations is of significant importance in microseismic investigations because this information greatly contributes to the overall success of hydraulic-fracturing monitoring programs. Full-wavefield time-reverse imaging (TRI) using one or more wave-equation imaging conditions offers an effective methodology for locating surface-recorded microseismic events. Although, to be most beneficial in microseismic monitoring programs, the TRI procedure requires using accurate subsurface models that account for elastic media effects. We have developed a novel microseismic (extended) PS energy imaging condition that explicitly incorporates the stiffness tensor and exhibits heightened sensitivity to isotropic elastic model perturbations compared with existing imaging conditions. Numerical experiments demonstrate the sensitivity of the microseismic TRI results to perturbations in P- and S-wave velocity models. Zero-lag and extended microseismic source images computed at selected subsurface locations yield useful information about 3D P- and S-wave velocity model accuracy. Thus, we assert that these image volumes potentially can serve as the input into microseismic elastic velocity model building algorithms.

Imaging Alpine crust using ambient noise wave-equation tomography

SUMMARYWe present an improved crustal Vs model and Moho depth map using ambient noise wave-equation tomography. The so-called ‘ambient noise wave-equation tomography’ is a method to invert seismic ambient noise phase dispersion data based on elastic waveform simulation, which accounts for 3-D and finite-frequency effects. We use cross-correlations of up to 4 yr of continuous vertical-component ambient seismic noise recordings from 304 high-quality broad-band stations in the Alpine region. We use model LSP_Eucrust1.0 obtained from traditional ambient noise tomography as initial model, and we iteratively improve the initial model by minimizing frequency-dependent phase traveltime differences between the observed and synthetic waveforms of Rayleigh waves in the period range 10–50 s. We obtain the final model after 15 iterations with ∼65 per cent total misfit reduction compared to the initial model. At crustal depth, the final model significantly enhances the amplitudes and adjusts the shapes of velocity anomalies. At Moho and upper-mantle depth, the final model corrects an obvious systematic velocity shift of the initial model. The resulting isovelocity Moho map confirms a Moho step along the external side of the external crystalline massifs of the northwestern Alps and reveals underplated gabbroic plutons in the lower most crust of the central and eastern Alps. Ambient noise wave-equation tomography turns out to be a useful tool to refine shear wave velocity models obtained by traditional ambient noise tomography based on ray theory.

Inline Delayed-shot Migration in Tilted Elliptical-cylindrical Coordinates

Riemannian wavefield extrapolation, a wave-equation imaging method for efficient one-way propagation of 3D wavefields, is extended to 3D seismic surveys. I formulate an inline delayed-shot migration procedure in 3D tilted elliptical-cylindrical (TEC) coordinate systems. When inline coordinate tilt angles are well matched to the inline source ray parameters, the TEC coordinate extension affords accurate propagation of both steep-dip and turning-wave components important for successfully imaging complex geologic structure. I show that wavefield extrapolation in TEC coordinates is no more complicated than propagation in elliptically anisotropic media, which is easily handled by current finitedifference practice. I apply this approach to a realistic 3D wide-azimuth synthetic derived from a field Gulf of Mexico data set. The resulting images demonstrate the imaging advantages made possible through 3D RWE implementations, including the improved imaging of steeply dipping salt flanks, at reduced computational cost.

Wave equation imaging and adjoint-state inversion for micro-seismic monitoring

Wave equation imaging and adjoint-state inversion for micro-seismic monitoring

Extended imaging conditions for passive seismic data

Seismic monitoring at injection sites (e.g., CO2 sequestration, hydraulic fracturing) has become an increasingly common tool amongst oil and gas producers. The information obtained from these data is often limited to seismic event properties (e.g., location, initiation time, moment tensor), the accuracy of which greatly depends on the assumed or estimated elastic velocity models. However, estimating accurate 3D velocity models from passive array data remains a challenging problem. Extended imaging conditions (eICs) for passive wave-equation imaging algorithms represent a key step towards generating – and verifying – elastic velocity models. By extending imaging conditions away from zero-lag in time and space we can better evaluate the focusing of a given event based on the principle that waves focus at zero lag only when the velocity models are correct. We demonstrate that given an elastic medium and multi-component recordings, we can propagate and correlate microseismic P- and S-wavefield modes to compute eICs for P- and S- velocity perturbations. We observe that the maximum correlation deviates from the zero-lag in time and space for a P/S cross-correlation imaging condition when using an incorrect P- and/or S-wave velocity, and thus there is sensitivity to velocity error not observable when using individual wavefield components.

Extended wave-equation imaging conditions for passive seismic data

The information obtained from seismic monitoring at injection sites is often limited to seismic event properties (e.g., location, origin time, moment tensor), the accuracy of which greatly depends on the assumed or estimated elastic velocity models. Current traveltime-based methods to calibrate velocity models rely on having distinct arrivals. In surface seismic monitoring, however, the signal-to-noise ratio of arrivals is often too poor for pick-based methods, and migration techniques are needed. Wave-equation migration can image microseismic sources without the need for picking or knowledge of the onset time by reconstructing the wavefield through the model space and applying an imaging condition. We have devised extended imaging conditions for passive seismic wave-equation imaging algorithms, and we have found that they represent a key step toward verifying and updating elastic velocity models. By extending imaging conditions away from zero lag in time and space, we can better evaluate the focusing of a given event based on a self-consistency principle that direct arrival waves focus at zero lag only when the velocity models were correct. We have determined that given an elastic medium and multicomponent recordings, we can propagate and correlate microseismic P- and S-wavefield modes to compute extended imaging conditions that are sensitive to P- and S-wave velocity perturbations. The extended image volume is robust to sparsely sampled data and high levels of noise and is sensitive to velocity model errors. The extended image volume can be used to determine relative P- and S-wave velocity updates to improve focusing.

Renormalized nonlinear sensitivity kernel and inverse thin-slab propagator in T-matrix formalism for wave-equation tomography

We derived the renormalized nonlinear sensitivity operator and the related inverse thin-slab propagator (ITSP) for nonlinear tomographic waveform inversion based on the theory of nonlinear partial derivative operator and its De Wolf approximation. The inverse propagator is based on a renormalization procedure to the forward and inverse transition matrix scattering series. The ITSP eliminates the divergence of the inverse Born series for strong perturbations by stepwise partial summation (renormalization). Numerical tests showed that the inverse Born T-series starts to diverge at moderate perturbation (20% for the given model of Gaussian ball with a radius of 5 wavelength), while the ITSP has no divergence problem for any strong perturbations (up to 100% perturbation for test model). In addition, the ITSP is a non-iterative, marching algorithm with only one sweep, and therefore very efficient in comparison with the iterative inversion based on the inverse-Born scattering series. This convergence and efficiency improvement has potential applications to the iterative procedure of waveform inversion.

An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging

Locating micro-seismic events is of utmost importance in seismic exploration, especially when searching for unconventional oil and gas resources. The arrival-time-difference approach is the dominant source location method currently used in the field of micro-seismic source location. However, micro-seismic events can be generated by any arbitrary rock movement and are often accompanied by interference noise. Recordings show characteristics of complicated wavelets and low signal-to-noise ratios. Under such conditions, conventional triangulation methods may have difficulty producing reliable locations; time-reversal imaging micro-seismic event location techniques are more promising. Locating micro-seismic events must be performed on-site for real-time monitoring of hydraulic fracturing. Introducing wave equation imaging techniques when locating micro-seismic events will increase the computation time, thus complicating real-time site monitoring. Therefore, the use of graphics processing unit (GPU) devices to accelerate time-reversal imaging micro-seismic event location technology becomes imperative. Three-dimensional synthetic data examples have demonstrated that the GPU-based time-reversal imaging micro-seismic event location method is typically 18 times faster than the central processing unit (CPU)-based implementation. The performance boost afforded by the GPU architecture allows us to locate micro-seismic events in 3D at a lower hardware cost and in less time than has been previously possible.

Inversion of seismic refraction and reflection data for building long-wavelength velocity models

Long-wavelength velocity model building is a nonlinear process. It has traditionally been achieved without appealing to wave-equation-based approaches for combined refracted and reflected waves. We developed a cascaded wave-equation tomography method in the data domain, taking advantage of the information contained in the reflected and refracted waves. The objective function was the traveltime residual that maximized the crosscorrelation function between real and synthetic data. To alleviate the nonlinearity of the inversion problem, refracted waves were initially used to provide vertical constraints on the velocity model, and reflected waves were then included to provide lateral constraints. The use of reflected waves required scale separation. We separated the long- and short-wavelength subsurface structures into velocity and density models, respectively. The velocity model update was restricted to long wavelengths during the wave-equation tomography, whereas the density model was used to absorb all the short-wavelength impedance contrasts. To improve the computation efficiency, the density model was converted into the zero-offset traveltime domain, where it was invariant to changes of the long-wavelength velocity model. After the wave-equation tomography has derived an optimized long-wavelength velocity model, full-waveform inversion was used to invert all the data to retrieve the short-wavelength velocity structures. We developed our method in two synthetic tests and then applied it to a marine field data set. We evaluated the results of the use of refracted and reflected waves, which was critical for accurately building the long-wavelength velocity model. We showed that our wave-equation tomography strategy was robust for the real data application.

Hybrid Kinematic-Dynamic Approach to Seismic Wave-Equation Modeling, Imaging, and Tomography

Estimation of the structure response to seismic motion is an important part of structural analysis related to mitigation of seismic risk caused by earthquakes. Many methods of computing structure response require knowledge of mechanical properties of the ground which could be derived from near-surface seismic studies. In this paper we address computationally efficient implementation of the wave-equation tomography. This method allows inverting first-arrival seismic waveforms for updating seismic velocity model which can be further used for estimating mechanical properties. We present computationally efficient hybrid kinematic-dynamic method for finite-difference (FD) modeling of the first-arrival seismic waveforms. At every time step the FD computations are performed only in a moving narrowband following the first-arrival wavefront. In terms of computations we get two advantages from this approach: computation speedup and memory savings when storing computed first-arrival waveforms (it is not necessary to make calculations or store the complete numerical grid). Proposed approach appears to be specifically useful for constructing the so-called sensitivity kernels widely used for tomographic velocity update from seismic data. We then apply the proposed approach for efficient implementation of the wave-equation tomography of the first-arrival seismic waveforms.

Validity of the Rytov Approximation in the Form of Finite-Frequency Sensitivity Kernels

The first-order (or linear) Rytov or Born approximation is the foundation for formulation of wave-equation tomography and waveform inversion, so the validity of the Rytov/Born approximation can substantially affect the applicability of these theories. However, discussions and research reported in literature on this topic are insufficient or limited. In this paper we introduce five variables in scattering theory to help us discuss conditions under which the Rytov approximation, in the form of the finite frequency sensitivity kernels (RFFSK), the basis of waveform inversion and tomography, is valid. The five variables are propagation length L, heterogeneity scale a, wavenumber k, anisotropy ratio ξ, and perturbation strength ɛ. Combined with theoretical analysis and numerical experiments, we conclude that varying the conditions used to establish the Rytov approximation can lead to uninterpretable or undesired results. This conclusion has two consequences. First, one cannot rigorously apply the linear Rytov approximation to all theoretical or practical cases without discussing its validity. Second, the nonlinear Rytov approximation is essential if the linear Rytov approximation is not valid. Different from previous literature, only phase (or travel time) terms for the whole wavefield are discussed. The time shifts of two specific events between the background and observed wavefields measured by cross-correlation will serve as a reference for evaluation of whether the time shifts predicted by the FFSKs are reasonably acceptable. Significantly, the reference “cross-correlation” should be regarded as reliable only if the condition “two specific similar signals” is satisfied. We cannot expect it to provide a reasonable result if this condition is not met. This paper reports its reliability and experimental limitations. Using cross-correlation (CC) samples as the X axis and sensitivity kernel (SK) or ray tracing (RT) samples as the Y axis, a chart of cross validation among the three is established. According to experimental analysis for a random and deterministic medium, the RFFSK works well for diffraction and geometric optics enclosed by the boundaries ka > 1 and (Λ 1, Φ < 2π), where Λ is a diffraction variable and Φ is a variable used to measure scattering strength. Both Λ and Φ depend on the five variables mentioned above. However, a large perturbation will limit the applicability of the Rytov approximation and narrow the range in which it is valid, because larger ɛ results in higher-amplitude fluctuations and larger back-scattering. Therefore, large ɛ is important in making SK/RT predictions discrepant from correct CC measurements. We also confirm that the RFFSK theory has wider range of application than ray theory. All the experiments point out the importance, both practical and theoretical, of respecting the conditions used to establish the Rytov approximation.

Resolution limits for wave equation imaging

Formulas are derived for the resolution limits of migration-data kernels associated with diving waves, primary reflections, diffractions, and multiple reflections. They are applicable to images formed by reverse time migration (RTM), least squares migration (LSM), and full waveform inversion (FWI), and suggest a multiscale approach to iterative FWI based on multiscale physics. That is, at the early stages of the inversion, events that only generate low-wavenumber resolution should be emphasized relative to the high-wavenumber resolution events. As the iterations proceed, the higher-resolution events should be emphasized. The formulas also suggest that inverting multiples can provide some low- and intermediate-wavenumber components of the velocity model not available in the primaries. Finally, diffractions can provide twice or better the resolution than specular reflections for comparable depths of the reflector and diffractor. The width of the diffraction–transmission wavepath is approximately λ at the diffractor location for the diffraction–transmission wavepath.

Integrated inversion using combined wave-equation tomography and full waveform inversion

Integrated inversion using combined wave-equation tomography and full waveform inversion

From tomography to full-waveform inversion with a single objective function

In full-waveform inversion (FWI), a gradient-based update of the velocity model requires an initial velocity that produces synthetic data that are within a half-cycle, everywhere, from the field data. Such initial velocity models are usually extracted from migration velocity analysis or traveltime tomography, among other means, and are not guaranteed to adhere to the FWI requirements for an initial velocity model. As such, we evaluated an objective function based on the misfit in the instantaneous traveltime between the observed and modeled data. This phase-based attribute of the wavefield, along with its phase unwrapping characteristics, provided a frequency-dependent traveltime function that was easy to use and quantify, especially compared to conventional phase representation. With a strong Laplace damping of the modeled, potentially low-frequency, data along the time axis, this attribute admitted a first-arrival traveltime that could be compared with picked ones from the observed data, such as in wave equation tomography (WET). As we relax the damping on the synthetic and observed data, the objective function measures the misfit in the phase, however unwrapped. It, thus, provided a single objective function for a natural transition from WET to FWI. A Marmousi example demonstrated the effectiveness of the approach.