Wavefield Information Research Articles

An algorithm for fast elastic wave simulation using a vectorized finite difference operator

Modern geophysical imaging techniques exploit the full wavefield information which can be simulated numerically. These numerical simulations are computationally expensive due to several factors, such as a large number of time steps and nodes, big size of the derivative stencil and huge model size. Besides these constraints, it is also important to reformulate the numerical derivative operator for improved efficiency. In this paper, we have introduced a vectorized derivative operator over the staggered grid with shifted coordinate systems. The operator increases the efficiency of simulation by exploiting the fact that each variable can be represented in the form of a matrix. This operator allows updating all nodes of a variable defined on the staggered grid, in a manner similar to the collocated grid scheme and thereby reducing the computational run-time considerably. Here we demonstrate an application of this operator to simulate the seismic wave propagation in elastic media (Marmousi model), by discretizing the equations on a staggered grid. We have compared the performance of this operator on three programming languages, which reveals that it can increase the execution speed by a factor of at least 2–3 times for FORTRAN and MATLAB; and nearly 100 times for Python. We have further carried out various tests in MATLAB to analyze the effect of model size and the number of time steps on total simulation run-time. We find that there is an additional, though small, computational overhead for each step and it depends on total number of time steps used in the simulation. A MATLAB code package, ’FDwave’, for the proposed simulation scheme is available upon request.

Sensitivity analysis of P-waves and S-waves to gas hydrate in the Shenhu area using OBS

Compared to towed streamers, ocean-bottom seismometers (OBS) obtain both S-wave data and richer wavefield information. In this paper, the induced polarization method is used to conduct wavefield separation on OBS data obtained from the Shenhu area in the South China Sea. A comparison of the changes in P- and S-waves, and a comprehensive analysis of geological factors within the area, enable analysis and description of the occurrence of natural gas hydrate in the study area. Results show an increase in P-wave velocity when natural gas hydrate exists in the formation, whereas the S-wave velocity remains almost constant, as S-waves can only propagate through the rock skeleton. Therefore, the bottom-simulating reflection (BSR) response of the P-wave is better than that of the S-wave in the frequency analysis profile. In a wide-angle section, the refractive wave of the hydrate layer is evident when using P-wave components but identification is difficult with S-wave components. This velocity model illustrates the sensitivity of P- and S-wave components to gas hydrate. The use of this polarization method and results of analysis provide technical and theoretical support for research on hydrate deposits and other geological features in the Shenhu area.

Solving the tensorial 3D acoustic wave equation: A mimetic finite-difference time-domain approach

Generating accurate numerical solutions of the acoustic wave equation (AWE) is a key computational kernel for many seismic imaging and inversion problems. Although finite-difference time-domain (FDTD) approaches for generating full-wavefield solutions are well-developed for Cartesian computational domains, several challenges remain when applying FDTD approaches to scenarios arguably best described by more generalized geometry. In particular, how best to generate accurate and stable FDTD solutions for scenarios involving grids conforming to complex topography or internal surfaces. We address these issues by developing a mimetic FDTD (MFDTD) approach that combines four key components: a tensorial 3D AWE, mimetic finite-difference (MFD) operators, fully staggered grids (FSGs), and MFD Robin boundary conditions (RBC). The tensorial formulation of the 3D AWE permits wave propagation to be specified on (semi-) analytically defined coordinate meshes designed to conform to complex domain boundaries. MFD operators allow for higher order FD stencils to be applied throughout the model domain, including the boundary region where implementing centered FD stencils can be problematic. The FSG approach combines wavefield information propagated on four complementary subgrids to ensure the existence of all wavefield gradients required for computing the tensorial Laplacian operator, and thereby avoids interpolation approximations. The RBCs are implemented with a flux-preserving mimetic boundary operator that forestalls introduction of nonphysical energy into the grid by enforcing underlying flux-conservation laws. After validating the 3D MFDTD scheme on a sheared Cartesian mesh, we generate 3D wavefield simulation examples for internal boundary (IB) and topographic coordinate systems. The numerical examples demonstrate that the MFDTD scheme is capable of providing accurate and low-dispersion impulse responses for scenarios involving distorted IB meshes conforming to water-bottom surfaces and topographic coordinate systems exhibiting 2.5 km of topographic relief and including steep (65°) slope angles.

The frequency-difference and frequency-sum acoustic-field autoproducts.

The frequency-difference and frequency-sum autoproducts are quadratic products of solutions of the Helmholtz equation at two different frequencies (ω+ and ω-), and may be constructed from the Fourier transform of any time-domain acoustic field. Interestingly, the autoproducts may carry wave-field information at the difference (ω+ - ω-) and sum (ω+ + ω-) frequencies even though these frequencies may not be present in the original acoustic field. This paper provides analytical and simulation results that justify and illustrate this possibility, and indicate its limitations. The analysis is based on the inhomogeneous Helmholtz equation and its solutions while the simulations are for a point source in a homogeneous half-space bounded by a perfectly reflecting surface. The analysis suggests that the autoproducts have a spatial phase structure similar to that of a true acoustic field at the difference and sum frequencies if the in-band acoustic field is a plane or spherical wave. For multi-ray-path environments, this phase structure similarity persists in portions of the autoproduct fields that are not suppressed by bandwidth averaging. Discrepancies between the bandwidth-averaged autoproducts and true out-of-band acoustic fields (with potentially modified boundary conditions) scale inversely with the product of the bandwidth and ray-path arrival time differences.

A new high-order finite volume method for 3D elastic wave simulation on unstructured meshes

In this paper, we proposed a new efficient high-order finite volume method for 3D elastic wave simulation on unstructured tetrahedral meshes. With the relative coarse tetrahedral meshes, we make subdivision in each tetrahedron to generate a stencil for the high-order polynomial reconstruction. The subdivision algorithm guarantees the number of subelements is greater than the degrees of freedom of a complete polynomial. We perform the reconstruction on this stencil by using cell-averaged quantities based on the hierarchical orthonormal basis functions. Unlike the traditional high-order finite volume method, our new method has a very local property like DG and can be written as an inner-split computational scheme which is beneficial to reducing computational amount. Moreover, the stencil in our method is easy to generate for all tetrahedrons especially in the three-dimensional case. The resulting reconstruction matrix is invertible and remains unchanged for all tetrahedrons and thus it can be pre-computed and stored before time evolution. These special advantages facilitate the parallelization and high-order computations. We show convergence results obtained with the proposed method up to fifth order accuracy in space. The high-order accuracy in time is obtained by the Runge–Kutta method. Comparisons between numerical and analytic solutions show the proposed method can provide accurate wavefield information. Numerical simulation for a realistic model with complex topography demonstrates the effectiveness and potential applications of our method. Though the method is proposed based on the 3D elastic wave equation, it can be extended to other linear hyperbolic system.

Applicability of SAR-based wave retrieval for wind–wave interaction analysis in the fetch-limited Baltic

ABSTRACTIn this article, a method for the detection of wave field parameters from synthetic aperture radar (SAR) imagery in the fetch-limited Baltic Sea is presented. Over the Baltic Sea region, common southwest (SW) and west (W) winds induce steep waves with shorter wavelengths compared with ocean waves. Thus, with the use of previous SAR sensors (e.g. ENVISAT/ASAR), it was not possible to detect individual waves and retrieve image wave number spectra. Since the year 2007, when TerraSAR-X (TS-X) reached its orbit, high spatial resolution data is available for measuring the sea-state parameters: the individual waves up to 30 m wavelength and their refraction can be distinguished. The main objective of this work was to demonstrate the capability of detecting wave field parameter from (TS-X) imagery in the Baltic Sea. The wave field parameters obtained from the SAR imagery were compared with in situ measurements and the Simulating WAves Nearshore (SWAN) wave model. The comparison of SAR-based wave field information with buoy measurements showed high agreement in case of wave propagation direction (r = 0.95) and wavelength (r = 0.83). A significant correlation is also seen between SWAN- and SAR-derived wave propagation direction (r = 0.87) and wavelengths (r = 0.91). With the case studies, it is shown that SAR data enables one to detect land shadow effects and small-scale wave field variations in the coastal zone. It was shown that SAR data is also valuable for improving and interpreting the wave model results. In consequence of common slanting fetch cases over the Baltic Sea region, it was demonstrated that the peak wave directions differ from the mean wind directions up to 43°.

Locating material defects via wavefield demixing with morphologically germane dictionaries

This article introduces a methodology for the detection and localization of structural defects in solid media using morphological demixing algorithms. The demixing algorithms are designed to separate spatiotemporal response data into two morphologically antithetical components: one contribution captures the spatially sparse and temporally persistent features of the medium’s response, while the other provides a representation of the dominant, globally smooth component as it would be observed in a defect-free medium. Within the demixing paradigm, we explore two methods: in the first, we cast the demixing task in terms of a group Lasso regularization problem with simply structured orthonormal dictionaries; the second method makes use of a more morphologically germane dictionary whose additional structure allows for the localization of defects whose signature may be highly elusive, for example, buried in noise or masked by competing features. After the demixing is complete, an automatic visualization tool highlights the regions associated with potential anomalies and outputs their local coordinates. Since the method does not invoke any knowledge of the material properties of the medium, or of its behavior in its pristine conditions, and is solely based on data processing of current wavefield information, it is endowed with significant model agnostic and baseline-free attributes. These properties are desirable in systems where there exists limited or unreliable a priori knowledge of the constitutive model, when the physical domain is highly heterogeneous or compromised by large damage zones, or when accurate baseline simulations are unavailable. The efficacy of the proposed method is evaluated against synthetically generated data and experimental data obtained using a scanning laser Doppler vibrometer.

Prismatic and full-waveform joint inversion

Prismatic wave is that it has three reflection paths and two reflection points, one of which is located at the reflection interface and the other is located at the steep dip angle reflection layer, so that contains a lot of the high and steep reflection interface information that primary cannot reach. Prismatic wave field information can be separated by applying Born approximation to traditional reverse time migration profile, and then the prismatic wave is used to update velocity to improve the inversion efficiency for the salt dame flanks and some other high and steep structure. Under the guidance of this idea, a prismatic waveform inversion method is proposed (abbreviated as PWI). PWI has a significant drawback that an iteration time of PWI is more than twice as that of FWI, meanwhile, the full wave field information cannot all be used, for this problem, we propose a joint inversion method to combine prismatic waveform inversion with full waveform inversion. In this method, FWI and PWI are applied alternately to invert the velocity. Model tests suggest that the joint inversion method is less dependence on the high and steep structure information in the initial model and improve high inversion efficiency and accuracy for the model with steep dip angle structure.

Verification of the helioseismic Fourier-Legendre analysis for meridional flow measurements

Measuring the Sun's internal meridional flow is one of the key issues of helioseismology. Using the Fourier-Legendre analysis is a technique for addressing this problem. We validate this technique with the help of artificial helioseismic data. The analysed data set was obtained by numerically simulating the effect of the meridional flow on the seismic wave field in the full volume of the Sun. In this way, a 51.2-hour long time series was generated. The resulting surface velocity field is then analyzed in various settings: Two $360^\circ \times 90^\circ$ halfspheres, two $120^\circ \times 60^\circ$ patches on the front and farside of the Sun (North and South, respectively) and two $120^\circ \times 60^\circ$ patches on the northern and southern frontside only. We compare two possible measurement setups: observations from Earth and from an additional spacecraft on the solar farside, and observations from Earth only, in which case the full information of the global solar oscillation wave field was available. We find that, with decreasing observing area, the accessible depth range decreases: the $360^\circ \times 90^\circ$ view allows us to probe the meridional flow almost to the bottom of the convection zone, while the $120^\circ \times 60^\circ$ view means only the outer layers can be probed. These results confirm the validity of the Fourier-Legendre analysis technique for helioseismology of the meridional flow. Furthermore these flows are of special interest for missions like Solar Orbiter that promises to complement standard helioseismic measurements from the solar nearside with farside observations.

Least-squares reverse-time migration with cost-effective computation and memory storage

Least-squares reverse-time migration (LSRTM), which involves several iterations of reverse-time migration (RTM) and Born modeling procedures, can provide subsurface images with better balanced amplitudes, higher resolution and fewer artifacts than standard migration. However, the same source wavefield is repetitively computed during the Born modeling and RTM procedures of different iterations. We developed a new LSRTM method with modified excitation-amplitude imaging conditions, where the source wavefield for RTM is forward propagated only once while the maximum amplitude and its excitation-time at each grid are stored. Then, the RTM procedure of different iterations only involves: (1) backward propagation of the residual between Born modeled and acquired data, and (2) implementation of the modified excitation-amplitude imaging condition by multiplying the maximum amplitude by the back propagated data residuals only at the grids that satisfy the imaging time at each time-step. For a complex model, 2 or 3 local peak-amplitudes and corresponding traveltimes should be confirmed and stored for all the grids so that multiarrival information of the source wavefield can be utilized for imaging. Numerical experiments on a three-layer and the Marmousi2 model demonstrate that the proposed LSRTM method saves huge computation and memory cost.

Full-azimuth anisotropic prestack time migration in the local-angle domain and its applications on fracture detection

Considering that geologic structures disturb prestack amplitude relationships, anisotropic migration is thus advocated not only for extracting azimuth-preserved common image gathers (CIGs), but also for preserving fracture-induced amplitude responses. However, most conventional anisotropic migration methods are hindered by their inefficiency in either modeling azimuthal traveltime variations at large offsets or characterizing subsurface reflections. Given that prestack time migration is widely applied for most practical purposes, we began with reformulations on a quartic traveltime formula, through which a new set of anisotropic parameters was developed. Then, an anisotropic migration method was established in the local-angle domain (LAD) for more reasonable uses of subsurface wavefield information. We also used a traveltime inversion scheme to estimate those anisotropic parameters required by anisotropic migration. Using this methodology on a physical model with a fracture medium, we derived better focused CIGs by thoroughly correcting the anisotropic effects of overburden. As a result, predicted properties of the fracture medium showed fewer interventions of geologic impacts. In a field example, a comprehensive study was performed on a deep carbonate reservoir to examine influences of different anisotropic migration algorithms on ultimate fracture prediction. Comparisons of the signal-to-noise ratio and agreements with formation microimage information reconfirmed the superiority of LAD anisotropic migration in recovering true properties of subsurface fractures, relative to routine methods (i.e., azimuth-sectored migration and anisotropic migration in the surface-offset domain).

Directional bending wave propagation in periodically perforated plates

We report on the investigation of wave propagation in a periodically perforated plate. A unit cell with double-C perforations is selected as a test article suitable to investigate two-dimensional dispersion characteristics, group velocities, and internal resonances. A numerical model, formulated using Mindlin plate elements, is developed to predict relevant wave characteristics such as dispersion, and group velocity variation as a function of frequency and direction of propagation. Experimental tests are conducted through a scanning laser vibrometer, which provides full wave field information. The analysis of time domain wave field images allows the assessment of plate dispersion, and the comparison with numerical predictions. The obtained results show the predictive ability of the considered numerical approach and illustrate how the considered plate configuration could be used as the basis for the design of phononic waveguides with directional and internal resonant characteristics.

Extracting depth information of 3-dimensional structures from a single-view X-ray Fourier-transform hologram.

We demonstrate how information about the three-dimensional structure of an object can be extracted from a single Fourier-transform X-ray hologram. In contrast to lens-based 3D imaging approaches that provide depth information of a specimen utilizing several images from different angles or via adjusting the focus to different depths, our method capitalizes on the use of the holographically encoded phase and amplitude information of the object's wavefield. It enables single-shot measurements of 3D objects at coherent X-ray sources. As the ratio of longitudinal resolution over transverse resolution scales proportional to the diameter of the reference beam aperture over the X-ray wavelength, we expect the approach to be particularly useful in the extreme ultraviolet and soft-X-ray regime.

Technique evolvement and improvement for better fracture characterization using prestack P-wave seismic data

Anisotropic migration has become a routine practice for fracture prediction because of its remarkable capability of suppressing effects of geologic structures while preserving anisotropy induced by fractures. The quality of input seismic gathers is of the first priority among other factors for fracture prediction. Azimuth-sectored migration loses considerable accuracy and reliability because of its coarse and inappropriate manner in dealing with seismic data. A full-azimuth anisotropic migration in surface offset domain provides better efficiency in using data information. Furthermore, full-azimuth anisotropic migration in local angle domain (subsurface angle domain) allows for accurate and reliable exploitation of seismic wavefield information. Finally, influences of the above three anisotropic migration methods on fracture prediction are prescribed with a field example in which validations with formation micro-imager (FMI) information demonstrate that full-azimuth anisotropic migration in local angle domain is the state-of-the-art technique for recovering true properties of subsurface fractures.

Gaussian beam prestack depth migration of converted wave in TI media

Abstract Increasing amounts of multi-component seismic data are being acquired on land and offshore because more complete seismic wavefield information is beneficial for structural imaging, fluid detection, and reservoir monitoring. S-waves are typically influenced more by anisotropy in a medium than are P-waves; as a result, the anisotropy cannot be ignored during the converted PS-wave imaging. Gaussian beam migration, an elegant and efficient depth migration method, is becoming a new topic in the study of PS-wave migration; its accuracy is comparable to that of wave-equation migration, and its flexibility is comparable to that of Kirchhoff migration. In this paper, we introduce an anisotropic Gaussian beam prestack depth migration (GB-PSDM) method for the converted PS-wave, in which the anisotropic media can be a transversely isotropic (TI) medium with a vertical or tilted symmetry axis. We present the PS-wave common shot gathers GB-PSDM imaging condition and derive the ray tracing of P- and SV-waves in two-dimensional TI media. The migration impulse responses of P- and SV-propagation modes in TI media with both vertical and tilted symmetry axes are presented. The results of numerical examples indicate that the method introduced here offers significant improvements in the quality of converted PS-wave imaging compared with an isotropic algorithm.

Elliptical moveout operator for data regularization in azimuthally anisotropic media

Data regularization by azimuthal moveout (AMO) is an important seismic processing step applied to minimize the deleterious effects of irregular and incomplete acquisition in complex geology. Using isotropic AMO operators on data acquired over azimuthally anisotropic media, though, can lead to poor regularization results due to mixing of wavefield information from neighboring traces with azimuthally varying velocity profiles. An elliptical moveout operator (EMO), representing an extension of isotropic AMO to elliptical azimuthally anisotropic media, is sensitive to variations in the magnitude and the orientation of velocity ellipticity. AMO and EMO operators can be applied to regularize data by moving wavefield information from traces acquired at neighboring offsets and midpoints to infill existing data holes. Unlike AMO, though, EMO operators also can be used in a data conditioning procedure to interpolate energy between seismic traces where input and output velocity profiles are azimuthally elliptical and isotropic, respectively. Resulting processed data volumes are approximately free of elliptical azimuthal anisotropy, as can be shown by comparing analytical traveltimes and numerically calculated wavefield arrivals. EMO thus represents a one-step regularization/conditioning procedure for elliptically azimuthally anisotropic media that is more consistent with wave-equation physics and yields more accurate results than when compared with those from isotropic processing and elliptical residual moveout operator static corrections.

Taming waveform inversion non-linearity through phase unwrapping of the model and objective functions

SUMMARY Traveltime inversion focuses on the geometrical features of the waveform (traveltimes), which is generally smooth, and thus, tends to provide averaged (smoothed) information of the model. On other hand, general waveform inversion uses additional elements of the wavefield including amplitudes to extract higher resolution information, but this comes at the cost of introducing non-linearity to the inversion operator, complicating the convergence process. We use unwrapped phase-based objective functions in waveform inversion as a link between the two general types of inversions in a domain in which such contributions to the inversion process can be easily identified and controlled. The instantaneous traveltime is a measure of the average traveltime of the energy in a trace as a function of frequency. It unwraps the phase of wavefields yielding far less non-linearity in the objective function than that experienced with conventional wavefields, yet it still holds most of the critical wavefield information in its frequency dependency. However, it suffers from non-linearity introduced by the model (or reflectivity), as reflections from independent events in our model interact with each other. Unwrapping the phase of such a model can mitigate this non-linearity as well. Specifically, a simple modification to the inverted domain (or model), can reduce the effect of the model-induced non-linearity and, thus, make the inversion more convergent. Simple numerical examples demonstrate these assertions.

Exploring the potentials and limitations of the time-reversal imaging of finite seismic sources

Abstract. The characterisation of seismic sources with time-reversed wave fields is developing into a standard technique that has already been successful in numerous applications. While the time-reversal imaging of effective point sources is now well-understood, little work has been done to extend this technique to the study of finite rupture processes. This is despite the pronounced non-uniqueness in classic finite source inversions. The need to better constrain the details of finite rupture processes motivates the series of synthetic and real-data time reversal experiments described in this paper. We address questions concerning the quality of focussing in the source area, the localisation of the fault plane, the estimation of the slip distribution and the source complexity up to which time-reversal imaging can be applied successfully. The frequency band for the synthetic experiments is chosen such that it is comparable to the band usually employed for finite source inversion. Contrary to our expectations, we find that time-reversal imaging is useful only for effective point sources, where it yields good estimates of both the source location and the origin time. In the case of finite sources, however, the time-reversed field does not provide meaningful characterisations of the fault location and the rupture process. This result cannot be improved sufficiently with the help of different imaging fields, realistic modifications of the receiver geometry or weights applied to the time-reversed sources. The reasons for this failure are manifold. They include the choice of the frequency band, the incomplete recording of wave field information at the surface, the excitation of large-amplitude surface waves that deteriorate the depth resolution, the absence of a sink that should absorb energy radiated during the later stages of the rupture process, the invisibility of small slip and the neglect of prior information concerning the fault geometry and the inherent smoothness of seismologically inferred Earth models that prevents the beneficial occurrence of strong multiple-scattering. The condensed conclusion of our study is that the limitations of time-reversal imaging – at least in the frequency band considered here – start where the seismic source stops being effectively point-localised.

Source localization with minimum variance distortionless response for spherical microphone arrays

To improve localization accuracy, the spherical microphone arrays are used to capture high-order wavefield information. For the far field sound sources, the array signal model is constructed based on plane wave decomposition. The spatial spectrum function is calculated by minimum variance distortionless response (MVDR) to scan the three-dimensional space. The peak values of the spectrum function correspond to the directions of multiple sound sources. A diagonal loading method is adopted to solve the ill-conditioned cross spectrum matrix of the received signals. The loading level depends on the alleviation of the ill-condition of the matrix and the accuracy of the inverse calculation. Compared with plane wave decomposition method, our proposed localization algorithm can acquire high spatial resolution and better estimation for multiple sound source directions, especially in low signal to noise ratio (SNR).

Direct measurement of the Kirkwood-Rihaczek distribution for the spatial properties of a coherent light beam

We present direct measurement of Kirkwood-Rihaczek (KR) distribution for spatial properties of coherent light beam in terms of position and momentum (angle) coordinates. We employ a two-local oscillator (LO) balanced heterodyne detection (BHD) to simultaneously extract distribution of transverse position and momentum of a light beam. The two-LO BHD could measure KR distribution for any complex wave field (including quantum mechanical wave function) without applying tomography methods (inverse Radon transformation). Transformation of KR distribution to Wigner, Glauber Sudarshan P- and Husimi or Q- distributions in spatial coordinates are illustrated through experimental data. The direct measurement of KR distribution could provide local information of wave field, which is suitable for studying particle properties of a quantum system. While Wigner function is suitable for studying wave properties such as interference, and hence provides nonlocal information of the wave field. The method developed here can be used for exploring spatial quantum state for quantum mapping and computing, optical phase space imaging for biomedical applications.