Quaternion SVD Methods for the Extraction of Rayleigh Waves

Separation of seismic wavefields is an important topic in geophysical research. One of the most challenging targets is the correct and complete separation of the body wavefield from the surface wavefield. When surface waves and body waves are superposed in the offset-time domain, the part of the signal constituting surface waves and the part of signal constituting body waves must be inferred from prior information or directly from the data. In general, real seismic waveforms are determined by a large number of factors and the correct wavefield separation is an ambitious target. Applications of such separation are manifold. The extracted surface waves can be used as input of the inverse problem of determining the S-waves velocities of the near surface ground layers. Vice versa, surface waves can be extracted in order to be removed from the seismic section. In reflection seismic, for example, only body waves constitute the signal and surface waves are considered noise. In particular, this thesis focuses on the extraction of Rayleigh waves, which are a particular type of surface waves with elliptical polarization, i.e., they induce ground particles to move on elliptic trajectories. Such waves require multicomponent sensors to be correctly recorded, as elliptic motion needs intrinsically at least bidirectional recording sensors. In this thesis we will use mainly only two components: the x and the z directions. Rayleigh waves, in spite of their multicomponent nature, are usually processed component wise, repeating the analysis for the x, y and z components. In order to introduce a more compact and effective multicomponent approach, I propose in this thesis a quaternion approach for the extraction of Rayleigh waves. In mathematics, quaternions are a number system that extends the complex numbers. They were first introduced by Sir W. R. Hamilton in 1843 [5], and they are commonly used for calculations involving three-dimensional rotations such as in three-dimensional computer graphics. Each quaternion is constituted by a scalar component and three different imaginary components. The imaginary components can be imagined as a vector part. I will use quaternions in order to represent the multicomponent signals recorded on triaxial geophones. In the geophysical field, quaternion algebra has been used to implement multicomponent Wiener deconvolution, to perform multicomponent velocity analysis via the implementation of quaternion matched filters, and to detect the wavefronts arrivals directions via the implementation of a biquaternion version of the MUSIC algorithm. The structure of the thesis is outlined in the following paragraphs. The first chapter is constituted by this introduction. In chapter two the theoretical propagation of the free modes of Rayleigh waves in elastic media is introduced. Two simple examples of ground are considered: the former is a homogeneous half space, which shows no dispersion, the latter is a fluid layer over a half space, which is one of the simplest ground model which shows dispersion characteristics; for this dispersive ground model, the dispersion curves equations are derived. In chapter three, first SVD is introduced, then a polarization method is described. This method uses a parameter inferred by the SVD analysis (the inverse ellipticity parameter) in order to determine the degree of ellipticity of the signal for local time windows. In order to improve the discrimina- tion between circularly/elliptically polarized signals and incoherent noise a filtering process that changes the phase-difference between two spatial components of the data is applied to the polarization maps. Furthermore, in this chapter the two datasets that are used in this thesis are introduced. They are a synthetic dataset generated for a 1D elastic ground model, and a real dataset acquired in China. In chapter four, the dispersion characteristics are derived for the synthetic and the real dataset. The dispersion characteristics are described by two dispersion curves, respectively the phase velocity dispersion curve and the group velocity dispersion curve. Phase velocity is the velocity of the wave ripples, and group velocity is the velocity of the wave packet. The phase velocity dispersion curves are determined in the frequency-(ray parameter) domain, also called f-p domain. No equivalent domain exists in the standard processing for computing the group velocity dispersion curves. I implemented an ad hoc transform in order to determine the group velocity and I called it frequency-(group slowness) transform. In this chapter an effective method for separating a mode of the Rayleigh wave will also be shown. It consists in the muting in the f-p domain of the area outside of a selected zone surrounding the mode. The computed dispersion curves are compared with the theoretical dispersion curves for the synthetic data. For the real data two modes of Rayleigh waves are identified. In chapter five quaternions and matrices of quaternions are introduced. The main definitions, properties, and different representations of quaternions are introduced. In addition, the unit quaternion representation of rotation and its relation with the 3 x 3 standard rotation matrix are described. Regarding quaternion matrices, the complex adjoint, and the definition of rank and SVD are introduced. In chapter six the representation of Rayleigh waves with quaternion SVD is studied. Firstly an example of signal of quaternion rank 1 is introduced, its form and its polarization characteristics are described. Next three test signals with different polarizations are introduced. They represent the three cases of circular, elliptic and linear polarizations. The different representations in quaternion autoimages for the three polarizations are examined. In particular, the first quaternion autoimage truncation of the elliptic test signal is shown to have circular polarization. I called this phenomenon circularization of the elliptic wave. Next, I will briefly point out dispersion effects and I will explain the motive for the introduction of the filter bank in the first QSVD extraction method. In chapter seven, the extraction method using a filter bank and QSVD is described. The filter bank is used to separate different frequency bands of the signal. In this chapter, firstly FIR filters and filter banks are introduced; then, the specific filter bank used for the extraction method is described and the procedure of the extraction method is outlined; finally the results applied on the synthetic and the real datasets are exhibited. In chapter eight, I will present the extraction method that performs the QSVD analysis on local portions of the seismogram in the offset-time domain. The results of this method applied on the synthetic and the real datasets are exhibited. The final chapter will reveal conclusions. In this chapter the results of the different extraction methods will be evaluated, and a global view of the thesis is given

Genetic algorithm full-waveform inversion: uncertainty estimation and validation of the results

We cast the genetic algorithm-full waveform inversion (GA-FWI) in a probabilistic framework that through a multi-step procedure, allows us to estimate the posterior probability distribution (PPD) in model space. Since GA is not a Markov chain Monte Carlo method, it is necessary to refine the PPD estimated by GA (GA PPD) via a resampling of the model space with a Gibbs sampler (GS), thus obtaining the GA+GS PPDs. We apply this procedure to two acoustic 2D models, an inclusion model and the Marmousi model, and we find a good agreement between the derived PPDs and the varying resolution due to changes in the seismic illumination. Finally, we randomly extract several models from the so derived PPDs to start many local full-waveform inversions (LFWIs), which produce final high-resolution models. This set of models is then used to numerically estimate the final uncertainty (GA+GS+LFWI PPD). The multimodal and wide PPDs derived from the GA optimization, become unimodal and narrower after LFWI and, in the well illuminated parts of the subsurface, the final GA+GS+LFWI PPDs contain the true model parameters. This confirms the ability of the GA optimization in finding a velocity model suitable as input to LFWI

Estimation of acoustic macro models using a genetic full-waveform inversion: Applications to the Marmousi model

We present a stochastic full-waveform inversion that uses genetic algorithms (GA FWI) to estimate acoustic macro-models of the P-wave velocity field. Stochastic methods such as GA severely suffer the curse of dimensionality, meaning that they require unaffordable computer resources for inverse problems with many unknowns and expensive forward modeling. To mitigate this issue, we propose a two-grid technique, that is, a coarse grid to represent the subsurface for the GA inversion and a finer grid for the forward modeling. We applied this procedure to invert synthetic acoustic data of the Marmousi model. We show three different tests. The first two tests use as prior information a velocity model derived from standard stacking velocity analysis and differ only for the parameterization of the coarse grid. Their comparison shows that a smart parameterization of the coarse grid may significantly improve the final result. The third test uses a linearly increasing 1D velocity model as prior information, a layer-stripping procedure, and a large number of model evaluations. All the three tests return velocity models that fairly reproduce the long-wavelength structures of the Marmousi. First-break cycle skipping related to the seismograms of the final GA-FWI models is significantly reduced compared to the one computed on the models used as prior information. Descent-based FWIs starting from final GA-FWI models yield velocity models with low and comparable model misfits and with an improved reconstruction of the structural details. The quality of the models obtained by GA FWI + descent-based FWI is benchmarked against the models obtained by descent-based FWI started from a smoothed version of the Marmousi and started directly from the prior information models. The results are promising and demonstrate the ability of the two-grid GA FWI to yield velocity models suitable as input to descent-based FWI

Estimating velocity macro-models using stochastic full-waveform inversion

During my Ph.D. program, I have investigated two different topics. The first topic (major topic) addresses the issue of determining a suitable starting model for Full-Waveform Inversion (FWI); this topic has occupied the majority of my Ph.D. work. The second topic concerns determining a method to remove peg-legs from marine data sets acquired with a towed dual-sensor streamer. I worked on this subject during the first year of my Ph.D. In the following paragraphs, I outline my work on both of these subjects. Estimating velocity macro-models using stochastic full-waveform inversion I have developed a procedure that estimates an acoustic 2D macro-model of the subsurface using a genetic algorithm and that uses the information of the entire seismic wavefield in the objective function. The aim of this work is to demonstrate that such an estimated 2D macro-model is well-suited to act as the starting model for high-resolution gradient-based full-waveform inversion. High-resolution gradient-based full-waveform inversion (which is usually referred to simply as full-waveform inversion (FWI)) is an iterative local optimization method that exploits all information from the seismogram to produce a high-resolution image of the subsurface. In the last twenty years, FWI has received increasing interest from both the oil and gas exploration industry and academia. In spite of being a very promising method, FWI is limited by its local nature, i.e., it terminates in the nearest local minimum. For this reason, starting the inversion from a good first-guess model is a crucial factor. To the best of our knowledge, an efficient method to determine a reliable starting model for FWI has not yet been found. Industry and academia have developed a number of procedures that might provide an adequate starting model, but they are usually very time-consuming. Furthermore, the majority of these procedures requires the picking for the arrival time of a number of seismic events and, generally, picking is a time-consuming task, which is tedious, subject to interpretation, and prone to errors. The approach that I propose does not require a tedious picking procedure, and, at the same time, it is resistant to falling into local minima. It is well-known that a starting model for FWI is a model that lies in the valley of the global minimum of a certain misfit surface. Consequently, the basic idea behind our approach is to attack the local nature of FWI by developing a global optimization method. A global method is able to escape from a local minimum because it is not driven by local derivatives of a misfit functional. A strong limitation of stochastic methods in geophysical inversion problems, and especially in FWI, is the so-called curse of dimensionality. To mitigate this issue, I have developed a simple strategy to reduce the number of unknowns of the model space in the synthetic inversions: each model of such a model space contains only the low wave-numbers and, hence, can be referred to as a macro-model. Another issue faced in this work is the high computational-cost of the stochastic full-waveform inversion. To reduce the overall computational-cost, I limited the seismic propagation from 3D to 2D and I performed several tests on an area of the subsurface smaller in size than those typically used in FWI tests. The time-step and space-step may be changed to further reduce the computational cost of each single forward-modeling. The most popular stochastic methods in Exploration Geophysics are genetic algorithms (GAs) and the simulated annealing (SA) method. Another popular method that originated in global seismology, is the neighborhood algorithm (NA). In this work, I compared a specific implementation of GAs, the Adaptive Simulated Annealing, and NA using two analytic objective functions and using a 1D elastic FWI problem. GA resulted the best performing method among the three for high-dimensional model spaces (>40). Consequently, I selected the GA, and I employed this method for the 2D acoustic full-waveform inversions on synthetic seismic data. The synthetic tests have been performed on both a portion and the entire Marmousi model. The Marmousi model is a synthetic model characterized by an intense layering, several faults, folds, and velocity inversions. Because of its complexity, it has been used widely as benchmark to test FWI algorithms. The outcomes of the synthetic tests have been employed as starting models for local FWI. I proved the validity of my methodology by comparing the final outcome after local FWI, with a reference workflow started from a smoothed version of the Marmousi model. A two-step depeg-leg method for marine acquisitions with towed dual-sensor streamers I present a two-step method that predicts (prediction step) and attenuates (subtraction step) peg-leg reflections in pre-stack seismic data acquired with towed dual-sensor streamers. The towed dual-sensor streamer permits to separate the wavefield in the down-going and up-going components. The advantage of having both the up-going and down-going wavefields separately available is exploited in the method. In this method, a key role is played by the shaping deconvolution in local windows of the data, which is applied to the predicted peg-leg wavefield prior to the subtraction step. I show that for windows in which no peg-leg signal is present, the shaping deconvolution may alter the primary reflections and it must not by applied. Hence, I added an automatic control to the local shaping deconvolution that preserves the primary signals during the peg-leg removal. Such a procedure has been applied to a synthetic data set generated by the reflectivity method. Using this data set, I verified the validity of the method and that the control added to the procedure dramatically improves the final result

Separation of a single mode of Rayleigh waves using quaternion SVD in vertically-heterogeneous media

I propose a wave-separation method that separates a single Rayleigh-wave mode from body waves and from other modes, using quaternions to represent the multi-component seismic data recorded by an array of vector-sensors. This method decomposes the signal into narrow-frequency bands, which are subjected to both a velocity correction and a polarisation correction. The aim of these corrections is to reduce the mode of interest to a quasi-monochromatic wave packet with infinite apparent velocity and quasi-circular polarisation. Once written in quaternion notation, I refer to this wave packet as “quaternion brick”. It can be proved that this quaternion brick maps into the first quaternion eigenimage of the quaternion SVD (QSVD). I apply this method to vertically-heterogeneous elastic models. The method seems to correctly extract the selected Rayleigh-wave mode and to separate it from body waves. The separation between different modes is more challenging. It results that, when two modes interfere in the near-offset portion, a more accurate separation can be obtained by characterising the mode of interest in the far-offset portion and then extrapolating it in the near-offset portion (where the two modes overlap) using the properties of the right singular vector of the QSVD

Strain Accumulation Mechanisms in Unconsolidated Sediments during Compression

Pride and Berryman (2009) proposed a model to predict pressure dependence of effective elastic bulk modulus for unconsolidated sediments, by progressively allowing the creation of new contacts during compression. They assume that the gaps around rattlers are distributed according to a power law with distance, in addition, the model allows two different strain-accumulation mechanisms: linear or quadratic, the latter being associated with grain rotation (Goddard, 1990). We have observed that the model of Pride and Berryman can be simplified without losing its generality, assuming a flat distribution of gaps around rattlers, given appropriate values for the maximum gap. We have used this simplified model to study how the strain-accumulation mechanism affects the coordination number during isotropic compression. We tested our model on sand data from Zimmer (2003). We observed that the majority of the experimental trends lay between the pure linear and the pure quadratic accumulation trends. We conclude that the strain accumulation in unconsolidated sediments can be well described as a combination of the two mechanisms. We noted also that rotation affects larger grains (diameter approx. 500 micron) more than smaller grains (diameter approx. 100 micron)

Characterisation and extraction of a Rayleigh-wave mode in vertically heterogeneous media using quaternion SVD

We propose a method that identifies a mode of Rayleigh waves and separates it from body waves and from other modes, using quaternions to represent multi-component data. Being well known the abilities of quaternions to handle rotations in space, we use previous results derived from Le Bihan and Mars (2004) to prove that a Rayleigh-wave mode recorded by an array of vector-sensors can be approximated by a sum of trace-by-trace rotating time signals. Our method decomposes the signal into narrow-frequency bands, which undergo both a velocity correction and a polarisation correction. The aim of these corrections is to reduce the mode of interest to a quasi-monochromatic wave packet with infinite apparent velocity and quasi-circular polarisation. Once written in quaternion notation, we refer to this wave packet as quaternion brick. Based on theoretical considerations, we prove that this quaternion brick maps into the first quaternion eigenimage of the quaternion singular value decomposition. We apply this method to synthetic datasets derived from two vertically heterogeneous models to extract the fundamental mode and we prove that it is correctly separated from either a higher mode of propagation or body waves with negligible residual. Results are presented in both time-offset and frequency-phase slowness domains

Including plastic behaviour in the Preisach-Mayergoyz space to find static and dynamic bulk moduli in granular media

We propose a modification of the Preisach-Mayergoyz (PM) space model to account for plastic deformation of a heterogeneous medium subjected to hysteretic non-linear elasticity. The PM space models the heterogeneous medium as a set of hysteretic nonlinear mesoscopic units which behave as switches that expand (turn on) and contract (turn off) at different pressures. The density distribution of these units describes the elastic behaviour of the medium. The PM model accounts for hysteresis but not for plastic deformation. We modify the model to include the plastic deformation by allowing the units to expand (turn on) at negative pressures. We implemented the elasto-plastic PM model using a discretized representation according to (Guyer et al. 1997). We tested this model on two loading cycles of a Gulf of Mexico beach-sand sample (Zimmer 2003). We compare the classical PM to our elasto-plastic PM and we highlight the increased ability to predict the dynamic bulk modulus as well as the enhancement in the representation of the effective stress-strain path, while being consistent with the original PM model

Comparisons of recent global optimization algorithms: tests on analytic objective functions and residual statics corrections

We compare the performance of six recent global optimization algorithms: Imperialist Competitive Algorithm (ICA), Firefly Algorithm (FA), Water Cycle Algorithm (WCA), Whale Swarm Optimization (WSO), Fireworks Algorithm (FWA) and Quantum Particle Swarm Optimization (QPSO). These methods have been introduced in the last few years and have found very limited or no applications in geophysical exploration problems thus far. The methods are first tested on two multi-minima analytic objective functions often used to test optimization algorithms: The Rastrigin and the Schwefel functions. Then, they are compared on the residual statics corrections, which is a highly non-linear geophysical optimization problem. In particular, we are interested in testing the convergence capabilities of these methods as the number of unknown model parameters increases. The different approaches are compared against a standard implementation of the Particle Swarm optimization (PSO), that is a popular global search method. The tests on the analytical functions and on the residual statics corrections demonstrate that FA, WCA, WSO and FWA outperform the other approaches in solving multi-minima and high-dimensional optimization problems. Conversely, PSO and ICA show limited exploration capabilities and lower convergence rates with respect to the other approaches