Earth:Surface wave inversion

Seismic inversion involves the set of methods which seismologists use to infer properties through physical measurements. Surface-wave inversion is the method by which elastic properties, density, and thickness of layers in the subsurface are obtained through analysis of surface-wave dispersion. The entire inversion process requires the gathering of seismic data, the creation of dispersion curves, and finally the inference of subsurface properties.

Figure 1. Rayleigh vs Love Waves. The small arrows show particle motion. Particle motion in Love waves is parallel to the surface and normal to the direction of propagation. Displacement in Rayleigh waves occurs in a retrograde elliptical motion normal to the surface and parallel to the direction of wave propagation.

Surface waves

Surface waves are seismic waves that travel at the surface of the earth, along the air/earth boundary. Surface waves are slower than P-waves(compressional waves) and S-waves(transverse waves). Surface waves are classified into two basic types, Rayleigh waves and Love waves. Rayleigh waves travel in a longitudinal manner (the wave motion is parallel to the direction of wave propagation) with particle motion in a retrograde elliptical motion (Figure 1). The Rayleigh waves result from the interaction between P-waves and vertically polarized S-waves. Conversely, Love waves travel in a traverse manner (Figure 1) (the wave motion is perpendicular to the direction of wave propagation), consisting of horizontally polarized S-waves. In seismology, surface waves are collected along with other seismic data, but are traditionally considered noise and an impedance in interpreting deeper reflection and refraction information. Seismologists usually modify seismic equipment and experimental procedures to remove surface wave information from the data. Earthquake seismologists however require the information seismic surface waves provide and thus design their equipment to amplify and gather as much information on these waves as possible. The work by early earthquake seismologists to extract substantial information from surface wave data was the basis for surface wave inversion theory.

Figure 2. Wavelength vs depth. Longer wavelength penetrates deeper.

Dispersion

The usefulness of surface waves in determining subsurface elastic properties arises from the way in which they disperse. Dispersion (geology) is the way in which surface waves spread out as they travel across the surface of the earth. Basically, if ten waves travel along the surface of the earth at the same speed, there is no dispersion. If several of the waves start to travel faster than the others, dispersion is occurring. Surface waves of varying wavelengths penetrate to different depths (Figure 2) and travel at the velocity of the mediums they are travelling through. Figure 2 was generated by plotting the amplitude of surface waves against depth. This was done for two different wavelengths. Both waves have the same total energy, but the longer wavelength has its energy spread out over a larger interval. If earth materials’ elastic parameters yield higher velocities with depth, longer wavelength surface waves will travel faster than those with shorter wavelengths. The variation of velocities with wavelength makes it possible to infer critical information about the subsurface. Dobrin (1951) uses a water disturbance example to illustrate the phenomenon that longer wavelengths tend to travel faster. This increase in speed with wavelength is seen for both group velocities and phase velocities. A wave group consists of waves at varying wavelengths and frequencies. Individual waves of a wave group are usually generated at the same time, but tend to spread out within the group because each wavelet travels at a different speed. A group velocity is basically the speed at which a wave group travels. A phase velocity is the speed at which an individual wave travels, having its own characteristic wavelength and frequency. Fourier theory tells us that a sharp impulse is made up of infinite frequency content in phase at one point. If each frequency travels at the same speed, that peak will remain intact. If each frequency travels at a different speed, that peak will spread out (Figure 3). This spreading out is dispersion. Phase and group velocity are both dependent on wavelength and are related by the equation

[math]\displaystyle{ V_{\mathrm{group}}=V_{\mathrm{phase}}-\lambda \frac{\delta V_{\mathrm{phase}}}{\delta \lambda} }[/math]

where V_group is the group velocity, V_phase is the phase velocity, and λ is the wavelength. When attempting surface wave inversion, phase velocities are used more often than group velocities because it is easier to create a dispersion curve of phase velocities. A dispersion curve is a plot of velocity versus frequency or wavelength. After the dispersion curve has been generated, a surface wave inversion process is performed to calculate the subsurface elastic properties. The accuracy of the dispersion curve is crucial in obtaining the correct subsurface elastic parameters from inversion.

Figure 3. Wavelengths of different frequencies spread out over time.

Elastic Properties

Elastic properties of the earth are those properties which affect the propagation of elastic waves. These properties are Lamé parameters and are used to relate stress to strain in isotropic media through Hooke’s law. Density is also related to elastic parameters through velocity equations for compressional and shear waves.

Data Gathering

Two main data gathering techniques are employed in gathering surface wave information. The two methods are spectral analysis of surface waves (SASW) and multi-channel analysis of surface waves (MASW). These techniques use either passive or active sources. Passive sources are simply ambient noise, while active sources include traditional seismic sources such as an explosive device or a steel plate being hit with a hammer. Overall, passive energy sources usually require more time when data gathering than active energy. Ambient noise is also more useful when it comes from random directions. The spectral analysis surface wave (SASW) technique requires the use of a spectral analyzer and at least two geophones. The spectral analyzer is used to study the frequency and phase of signals being recorded by the geophones. An expanding spread array is useful in minimizing the near field effects of surface waves. An increase in offset distance will result in more time for the waves to reach each geophone, giving the longer wavelengths more time to disperse. The shot gather is modified to minimize the influence of body waves. As the data is gathered, the spectral analyzer is able to generate the dispersion curves for the survey area in real time. The multi-channel analysis of surface waves (MASW) technique can be performed similar to a traditional seismic acquisition whereby there is a geophone spread that is acquiring seismic data. The resulting data is processed by picking out the surface wave arrivals from the acquired distance vs. time plot. Based on the distance vs. time plot, the dispersion curve is created.

Dispersion curves

The process of creating dispersion curves from raw surface wave data (distance vs. time plot) can be performed using five transformation processes. The first is known as the wave-field transformation (τ-p transformation), first performed by McMechan and Yedlin (1981). The second is a 2-dimensional wave-field transform (f-k transformation) performed by Yilmaz (1987). The third is a wave-field transform base on phase shift, performed by Park et al. (1998). The fourth is a modified wave-field transform base on frequency decomposition and slant stacking, performed by Xia et al. (2007). The fifth is a high-resolution Linear Radon transformation performed by Luo et al. (2008). In performing a wave-field transformation, a slant stack is done, followed by a Fourier transform. The way in which a Fourier transform changes x-t data into x-ω (ω is angular frequency) data shows why phase velocity dominates surface wave inversion theory. Phase velocity is the velocity of each wave with a given frequency. The modified wavefield transform is executed by doing a Fourier transform first before a slant stack. Slant stacking is a process by which x-t (where x is the offset distance, and t is the time) data is transformed into slowness versus time space. A linear move (similar to normal move out (NMO)) out is applied to the raw data. For each line on a seismic plot, there will be a move out that can be applied that will make that line horizontal. Distances are integrated for each slowness and time composition. This is known as a slant stack because each value for slowness represents a slant in x-t space and the integration stacks these values for each slowness.

Modified wavefield transform

A Fourier transform is applied to raw surface wave data plotted x-t. u(x,t) represents the entire shot gather, and the Fourier transformation results in U(x,ω).

[math]\displaystyle{ U(x,\omega)= \int u(x,t)e^{-i \omega t} \,dt }[/math]

U(x,ω) is then deconvolved and can be expressed in terms of phase and amplitude.

[math]\displaystyle{ U(x,\omega)= P(x,\omega)A(x,\omega) }[/math]

where P(x,ω) is the phase portion of the equation that holds information containing the waves’ dispersion properties, including arrival time information and A(x,ω) is the amplitude portion that contains data pertaining to the attenuation and spherical divergence properties of the wave. Spherical divergence is the idea that as a wave spreads out, the energy in the wave spreads out over the surface of the waveform. Since P(x,ω) contains the dispersion property information,

[math]\displaystyle{ U(x,\omega)= e^{-i \Phi x} A(x,\omega) }[/math]

where Φ=ω/c_ω, ω is the frequency in radians, and c_ω is the phase velocity for frequency ω. This data can then be transformed to give velocity as a function of frequency:

[math]\displaystyle{ V(\omega,\Phi)= \int e^{i \Phi x} \frac{U(x,\omega)}{|U(x,\omega)|}\,dx }[/math]

This will yield a dispersion curve showing a variety of frequencies travelling at different phase velocities.

The surface wave inversion process is the act of inferring elastic properties such as density, shear wave velocity profile, and thickness from dispersion curves created. There are many methods (algorithms) that have been utilized to perform inversion including:

• Multilayer dispersion computation
• Least squares curve fitting program
• Knopoff’s method
• Direct search algorithm
• High frequency Rayleigh wave inversion
• Refraction microtremor method

Figure 4. Example of a dispersion curve where velocity increases with depth. The blue area represents experimental data, while the red line represents an experimental curve fit to the data.

Multilayer dispersion computation

Haskell (1953) first performed the multilayer dispersion computation. Haskell’s work has been the basis for much of the current surface wave inversion theory. Since Rayleigh waves are composed of P and S-waves and Love waves are composed of only S waves, Haskell derived the elastic wave equations for both P and S-waves. These equations were modified to show Rayleigh wave motion. After assuming a free surface boundary where no stresses or strains cross, the Rayleigh wave equation is simplified. Inputting different values for layer thicknesses, densities, and elastic parameters in the form of P and S wave velocities into the equation will yield a dispersion curve. Parameters can be modified to fit the derived dispersion curve to actual data (Figure 4).

Least squares curve fitting program

Dorman and Ewing (1962) came up with an algorithm based on Haskell’s earlier work. Their method used an iterative technique that enabled the user to input parameters and the computer to find which exact parameters best fit the experimental data.

Knopoff’s method

Knopoff’s method also uses Haskell’s equations to perform the surface wave data inversion, but it simplifies the equations for the fastest computation. The increased speed is mostly accomplished in programming as well as the lack of complex numbers in the calculations. In this algorithm, approximate layer thicknesses, compressional and shear velocities, as well as density values must be input for the model.

Direct search algorithm

The direct search algorithm matches a data driven model to the synthetic dispersion curve (Wathelet et al., 2004). This algorithm creates a theoretical dispersion curve by guessing parameters such as shear wave velocity, compressional wave velocity, density, and thickness. After the theoretical curve is created, the computer then attempts to match this theoretical curve with the actual (experimental) dispersion curve. The values of the parameters are picked at random, with different permutations, and repeated continuously until matching curves are achieved. In some cases, while running the algorithm, different values of shear and compressional velocities, density, and thickness might produce the same dispersion curve. The algorithm calculates a value known as the misfit value as it generates each theoretical dispersion curve. The misfit value is simply a measure of how the generated model stacks up to a true solution. Misfit is given by,

[math]\displaystyle{ Misfit = \sqrt{\sum_{i-\sigma}\frac{(x_{di}-x_{ci})^2}{\sigma^2_i n_F}} }[/math]

where x_di is the velocity of data curve at frequency f_i, x_ci is the velocity of the calculated curve at frequency f_i, σ_i is the uncertainty of the frequency samples considered and n_F is the number of frequency samples considered. If no uncertainty is provided, σ_i is replaced by x_di.

High frequency Rayleigh wave inversion

The high frequency Rayleigh wave inversion performed by Xia et al. (1999) analyzed the earth using Knopoff’s method. By varying different properties used in creating the dispersion curve, it was discovered that different earth properties had significantly different effects on phase velocities. Changing the S-wave velocity input has a dramatic impact on Rayleigh wave phase velocities at high frequencies (greater than 5 Hz). A change in S-wave velocity of 25% changes the Rayleigh wave velocity by 39%. Conversely, P-wave velocity and density have a relatively small impact on Rayleigh wave phase velocity. A change in density of 25% will cause a less than 10% change in surface wave velocity. A change in P-wave velocity will have even less effect (3%).

Microtremor method

The final inversion method, the refraction microtremor (ReMi) technique, makes use of a computer algorithm that forward models normal mode dispersion data obtained from a survey. This method uses regular P-wave and simple refraction acquisition equipment, and does not require an active source, hence the name. Pullammanapellil et al. (2003) used this method to accurately match the S-wave profile of the ROSRINE borehole drilled. The ReMi method accurately matched the overall shear wave velocity profile, but cannot match the detail provided by the shear velocity well log. The discrepancy in overall detail should have no effect in evaluating the subsurface.

Advantages/Disadvantages of Surface Wave Inversion

There are many advantages to using surface waves to image the subsurface. For one, surface wave inversion readily images low-velocity zones. Refraction methods cannot see low-velocity zones because such a zone would bend the traversing wave deeper instead of towards the surface. Surface wave inversion is also non-invasive as well as cost effective. There are a few disadvantages to this method as well. The resolution of the surface wave inversion method is not nearly as resolved as a seismic collection done in a wellbore. There is also the possibility for non-unique solutions to dispersion curves (several sets of parameters can yield the same dispersion curve). Additionally, the presence of multiple modes may exist and leak into the targeted mode that is being inverted for.

Conclusion

Surface wave inversion is becoming a valuable tool in evaluating the near subsurface. Surface waves found on seismograms can now be a useful by product of seismic exploration surveys instead of a waste product. Furthermore, it is more budget friendly because the use of an active energy source is not needed. Also, it is useful in detecting low velocity zones in the subsurface that are undetectable by refraction methods. It is most effective in estimating shear velocity, density, and thickness of subsurface profiles.

See also

• Spectral-Analysis-of-Surface-Waves
• Full waveform inversion

Uncited References

Foti, S., Comina, C., Boiero, D., Socco, L. V., 2009, Non-uniqueness in surface-wave inversion and consequences on seismic site response analyses: Soil Dynamics and Earthquake Engineering, v. 29, p. 982-993.

Kennett, B.L.N., 1976, The inversion of surface wave data: Pure and Applied Geophysics, v. 114, p 747-751.

Luke, B., Calderon-Macias, C., 2007, Inversion of seismic surface wave data to resolve complex profiles: Journal of Geotechnical and Geoenvironmental Engineering, v. 133, p. 155-165.

Lai, C. G., Foti, S., and Rix, G. J., 2005, Propagation of data uncertainty in surface wave inversion: Journal of Environmental & Engineering Geophysics, v. 10, p. 219-228.

Park, C., Miller, R., Laflen, D., Neb, C., Ivanov, J., Bennet, B., Huggins, R., 2004, Imaging dispersion curves of passive surface waves: SEG Expanded Abstracts, v. 23.

Supranata, Y. E., Kalinski M. E., Ye, Q., 2007, Improving the uniqueness of surface wave inversion using multiple-mode dispersion data: International Journal of Geomechanics, v. 7, p. 333-343.

Xia, J., Miller, R.D., Yixian, X., Yinhe, L., Chao, C., Jiangping, L., Ivanov, J., Zeng, C., 2009, High Frequency Rayleigh-Wave Method: Journal of Earth Science, v. 20, p. 563-579.

Yamanaka, H., Ishida, H., (1996). Application of genetic algorithms to an inversion of surface dispersion data: Bulletin of the Seismological Society of America, v. 86, p. 436-444.

Kallivokas, L.F., Fathi, A., Kucukcoban, S., Stokoe II, K.H., Bielak, J., Ghattas, O., (2013). Site characterization using full waveform inversion: Soil Dynamics and Earthquake Engineering, v. 47, p. 62-82.

Foti, S., Lai, C.G., Rix, G.J., and Strobbia, C., (2014). Surface Wave Methods for Near-Surface Site Characterization, CRC Press, Boca Raton, Florida (USA), 487 pp., ISBN:9780415678766 <https://www.crcpress.com/product/isbn/9780415678766>