Microtremor H/V spectral ratio (MHVSR) has gained popularity to assess the dominant frequency of soil sites. It requires measurement of ground motion due to seismic ambient noise at a site and a relatively simple processing. Theory asserts that the ensemble average of the autocorrelation of motion components belonging to a diffuse field at a given receiver gives the directional energy densities (DEDs) which are proportional to the imaginary parts of the Green’s function components when both source and receiver are the same point and the directions of force and response coincide. Therefore, the MHVSR can be modeled as the square root of 2 × ImG11/ImG33, where ImG11 and ImG33 are the imaginary parts of Green’s functions at the load point for the horizontal (sub-index 1) and vertical (sub-index 3) components, respectively. This connection has physical implications that emerge from the duality DED force and allows understanding the behavior of the MHVSR. For a given model, the imaginary parts of the Green’s functions are integrals along a radial wavenumber. To deal with these integrals, we have used either the popular discrete wavenumber method or the Cauchy’s residue theorem at the poles that account for surface waves normal modes giving the contributions due to Rayleigh and Love waves. For the retrieval of the velocity structure, one can minimize the weighted differences between observations and calculated values using the strategy of an inversion scheme. In this research, we used simulated annealing but other optimization techniques can be used as well. This last approach allows computing separately the contributions of different wave types. An example is presented for the mouth of Andarax River at Almería, Spain.

Among the various procedures to obtain subsoil properties, the microtremor horizontal-to-vertical spectral ratio, MHVSR, proposed by Nakamura (1989), is quite popular because it is a simple and economic way to estimate the dominant frequency of local stratigraphy using seismic ambient vibrations. For a comprehensive review, see Bonnefoy-Claudet et al. (2006). The MHVSR was assumed transfer function. Consequently, the success was accompanied of some controversy regarding such interpretation. For instance, Arai and Tokimatsu (2004) assumed the excitation of Love and Rayleigh waves and assigned their participation amount. They proposed the use of average directional energy densities (DEDs) to compute MHVSR as the square root of a ratio of DEDs. The use of energies opened the door for theoretical approach. Albarello and Lunedei (2011) and Lunedei and Albarello (2015) computed the energy densities for surface waves from a random distribution of surface sources (DSS). On the other hand, the emergence of diffuse field theory (e.g., Weaver 1982, 1985; Lobkis and Weaver 2001; Campillo and Paul 2003; Shapiro and Campillo 2004) prompted the research of coda and noise properties in controlled experiments (e.g., Hennino et al. 2001). Margerin et al. (2009) used DFA to compute the reciprocal of MHVSR squared for a layered structure. They did not pointed out the relationship with Nakamura’s (1989) proposal of H/V, probably because the focus of paper was on energy partitions and on other theoretical issues.

Weaver (1982, 1985) established energy densities within a diffuse field in a homogeneous elastic full-space and described them as either to wave types or to degrees of freedom. This can be analytically extended to the layered half-space system. The DEDs are proportional to the imaginary part of Green’s function for given directions and correspond to the power injected into the system by a unit harmonic load. They are the sum of contributions from a radiation episode and from the waves returning to the source. This expresses a duality of the radiated energy due to forcing the system at a point, and the “projection” of the diffuse field energy at such point, by means of correlations that extract the appropriate quantities from a ubiquitous diffuse field.

We assume that the ambient seismic noise is diffuse. This assumption embodies a practical strategy: The MHVSR emerges from field measurements and is modeled for a given configuration of the medium. This last task is easier to say than to accomplish. Once the illumination is set, the multiple scattering and the inhomogeneity of the studied region control the characteristics of the diffuse field. The theoretical examples considered are a full-space, a half-space, a layered medium and a homogeneous inclusion. In what follows, a brief account is given.

Isotropy of illumination and equipartition of the energy allow for exact retrieval of Green’s functions (see Sánchez-Sesma and Campillo 2006). The isotropy implies uniform coverage of the illumination and is a necessary condition to get the Green’s functions even in the presence of heterogeneity (Sánchez-Sesma et al. 2006). The term “diffuse” applies to the theoretical settings that allow the exact retrieval of Green’s functions. It is usual to extend the concept (perhaps somewhat loosely) to fields engendered by multiple scattering and the interactions with inclusions and voids in which it is definitively a flux of energy typical of diffusion equations. On the other hand, equipartition means that the available energy in the phase space goes with equal average amounts, to all the possible states. Extending the results by Weaver (1982, 1985) Sánchez-Sesma et al. (2008) and Perton et al. (2009) pointed out that in addition to modes the equipartition allows description with degrees of freedom. These two aspects of equipartition lead to the interpretation of microtremor H/V spectral ratios or MHVSR in terms of the imaginary part of Green’s functions (Sánchez-Sesma et al. 2011b) using diffuse field assumption (DFA). Recent developments include applications of the theory to study H/V for earthquakes (Kawase et al. 2011) receivers at depth (Lontsi et al. 2015) and sites with lateral irregularity (Matsushima et al. 2014).

For layered medium, the Green’s functions are integrals in the horizontal wavenumber domain (Sánchez-Sesma et al. 2011b). The Cauchy’s residue theory leads to a fast computation and allows separation of the various wave components. Thus, the inversion of H/V to obtain the soil profile with only one station is possible with the DFA (see García-Jerez et al. 2016 and Piña-Flores et al. 2016). To mitigate non-uniqueness, joint inversion with dispersion curves gives good results.

Sánchez-Sesma and Campillo (2006) showed that it is possible to retrieve the Green’s function of dynamic elasticity in an elastic space averaging cross-correlations of a diffuse field. They constructed such an isotropic field with a uniform, set of plane waves coming from all directions with the same energy. The vector cases imply the existence of P and S waves. By counting modes, Weaver (1982) showed that in a homogeneous full-space the energy ratio is ξS/ξP = 2α3/β3 = 2R3 where ξS and ξP are the energy densities of S and P waves, and β and α are their propagation velocities.

The use of Somigliana’s representation theorem allowed Sánchez-Sesma et al. (2008) to confirm and generalize the previous result. They obtained the expression:

Equation 1 is of great generality; although it implies a homogeneous envelope, it is applicable to any type of heterogeneity, inclusions and voids (see Sánchez-Sesma et al. 2006). When source and receiver coincide, \({\mathbf{x}}_{\text{A}} = {\mathbf{x}}_{\text{B}} = {\mathbf{x}}\), the left-hand side of Eq. 1 is precisely the average of spectral density. This is a finite quantity. Assuming i = j = 1, the directional energy density in direction 1 is:

The right-hand side is also finite: The singularity of Green’s function is restricted to the real part; the imaginary part is smooth and represents the radiated energy into the medium by the unit harmonic load. To explore the connection with a classical solution, consider the Fourier transform of Stokes’ (1849) solution when the distance source–receiver is zero:

where ξi is the energy density in direction i. Clearly, ξ1 + ξ2 + ξ3 = ξ which shows that energy can be also described using degrees of freedom.

Weaver (1985) studied diffuse fields at the free surface of half-space and found the contributions of P, SV, SH and Rayleigh waves to surface energy densities. Sánchez-Sesma et al. (2008) and Perton et al. (2009) also included the free surface to have a half-space in 2D and 3D, respectively, and studied the transition between deep space and the free surface. In fact, the energy density per unit area associated with Rayleigh waves is:

Perton et al. (2009) computed the energy densities against depth in an elastic half-space. They considered the incidence of an equipartitioned set of elastic waves. Figure 1 shows the energy densities in terms of either states (due to incoming P, SV, SH and Rayleigh waves) or degrees of freedom. As one may suspect, the total energy density is the same. This extends the elementary results for the full-space to the half-space with free surface and firmly establishes the concept of directional energy density.

Fig. 1

Left Energy densities for the three orthogonal directions. Energies E1, E2, and E3 with black discontinuous line. Right Energy densities associated with incoming P, SV, SH, and Rayleigh waves. Both plots are against normalized depth x3/λR, where x3 = vertical coordinate pointing toward the depth, λR = wavelength of Rayleigh waves. Left gray continuous comes from the imaginary part of Green’s function components at the source from Weyl type integral representation. In the right panel, the cases of SV and Rayleigh waves appear with discontinuous lines (from Perton et al. 2009, with the permission of AIP Publishing). We normalize the amplitudes to match the energy density in deep space

These two aspects of equipartition are the key, the connection with deterministic results for partition of energy due to surface loads, and lead to the interpretation of microtremor H/V spectral ratios or MHVSR in terms of the imaginary part of Green’s functions (see Sánchez-Sesma et al. 2011b) using the diffuse field assumption (DFA).

For a normal load on the surface of an elastic half-space, Miller and Pursey (1955) computed the partition of the injected energy by the load among the various types of elastic waves. For a Poisson solid about 2/3 of the energy leaves the loaded point as Rayleigh surface waves. These results were serendipitously generalized by Weaver (1985) for all Poisson ratios, including the tangential load case. He was looking for energy densities at the free surface. Sánchez-Sesma et al. (2011a) pointed out that calculations based on the diffuse fields concepts allow obtaining the partition of the energy injected on the elastic half-space by a surface load. This interpretation allows drawn parallels between and DEDs and applied forces.

The microtremor horizontal-to-vertical spectral ratio, MHVSR, proposed by Nakamura (1989) is very popular because it is a simple, economic way to estimate the dominant frequency of a stratigraphy using seismic ambient vibrations. It was first interpreted as a transfer function, which is not, and therefore, its success was accompanied of some controversy regarding the interpretation of properties. It was suggested that H/V to be a proxy of ellipticity. Arai and Tokimatsu (2004) proposed the use of average directional energy densities (DEDs) to compute MHVSR as the square root of a ratio of DEDs. In the two calculations, the main peak of H/V is preserved. While Nakamura (1989) proposed the average of ratios, Arai and Tokimatsu (2004) used the ratio of averages. The energy densities provide the connection between field measurements and intrinsic properties of systems. The spectral ratio H/V using measured energy densities is:

which is an intrinsic property of the system (Sánchez-Sesma et al. 2011b). Numerical results were obtained for the surface of a layered medium (x = 0). An early calculation of the reciprocal of MHVSR squared for a layered structure using DFA was due to Margerin et al. (2009). When body waves are the dominant part of wavefield, Kawase et al. (2011) gave an approach closely related to Claerbout (1968). Considering the DFA in 3D, one can model the MHVSR for measurements at depth (Lontsi et al. 2015) and for sites with lateral heterogeneity (Matsushima et al. 2014).

For horizontally layered medium, the imaginary parts of Green’s functions imply integrals of the form:

where fSVV(k), fSVH(k), and fSH(k) are the Harkrider’s (1964) coefficients in terms of Haskell’s (1953) propagators. Here k = radial wavenumber. These coefficients result from the global matrix (Knopoff 1964) as well (see Sánchez-Sesma et al. 2011b). Direct computations require very small discretization step. To speed-up calculations, the Bouchon’s (2003) discrete wavenumber method (DWN) is an option. The use of Harvey’s (1981) locked mode approximation is providing accurate and faster calculations (Wu et al. 2017, Personal communication). As we want also to isolate the contributions of various waves, contour integration (García-Jerez et al. 2013) is used.

In Fig. 2, the results of Im[G11(0, 0; ω)] and Im[G33(0, 0; ω)] for a layer over a rigid half-space are shown. A layer has thickness of 125 m and shear wave velocity of 500 m/s. The half-space velocities are ten times those of layer (Poisson ratio is 0.25) and equal densities. The jumps of ImG11 reveal the efficient injection of body shear waves by a horizontal load. It is clear that the shear wave resonances mark the emergence of higher modes of 1D response and of Love and Rayleigh waves. The behavior of ImG33 has also deep physical meaning. For low-frequency spectrum, it shows a cutoff, but for larger frequencies, it behaves as a half-space with the properties of layer. For very short wavelengths, the layer is like a “mini half-space” with significant emission of Rayleigh waves. The H/V computed using Eq. 9 appears with black line in Fig. 2 while the result using generalized eigenfunctions (Margerin et al. 2009) is also shown. The agreement is very good.

Fig. 2

(Left) Imaginary part of Green’s functions G11 and G33 for a layer over a half-space. Black line corresponds to the horizontal response while dotted line represents the vertical one. (Right) H/V computed with Eq. 9. Dotted line was computed with Margerin’s (2009) method

After the initial tests, Dr. T. Yokoi of Building Research Institute of Tsukuba, Japan, proposed a blind test with a model studied at SESAME project (Bonnefoy-Claudet et al. 2004; Bard 2008). Figure 3 shows model N104, which has two layers with sharp variations. The excitation was synthetic noise from random sources. Figure 3 depicts H/V at the surface of model for synthetic noise plotted along with their deviations (Bonnefoy-Claudet et al. 2004). Theoretical results are also given.

Fig. 3

(Left) Stratigraphy profiles of model N104 studied within the SESAME Project. (Right) Black line shows H/V computed at the surface of model N104 excited with synthetic noise. The bars indicate the standard deviation. With black continuous line and gray dotted line, the results of the computations using the DFA with no damping (1/Q = 0) and with 1/Q = 0.005 (Q = 200), respectively

The effects of attenuation are present in real data, but using damping within the diffuse field assumption is not considered. Snieder et al. (2007) showed that Green’s function emerges from the response to random forcing for a variety of conditions, including the extreme case of a model fulfilling the diffusion equation. In practice, small amounts of damping allow stable results and the retrieved Green’s functions are consistent with observations (e.g., Lawrence and Prieto 2011).

The assumed model was perfectly elastic, and no damping was required in the modeling with DFA. The simulations were performed with a numerical method and assumed it considers some small attenuation. We did not pursue a detailed explanation but consider some scenarios. In fact, calculations were performed with no damping at all (1/Q = 0) and with damping h = 0.25% (1/Q = 0.005). This last value gives results that matched data best. We considered the blind test satisfactory. With no damping, the computed first peak is 1.23 times the simulated one.

For horizontal layering, a fast calculation of the Green’s functions is possible and it is well suited for inversion. For irregular layers with smooth lateral variations, in terms of the wavelength, the 1D model could be useful. Forward problem is nonlinear and depends on several uncorrelated parameters. Moreover, several sets of parameters may give the same H/V. Figure 4 shows two velocity profiles. They give rise to the same H/V curve, demonstrating clearly the non-unicity of the H/V.

Because of the non-uniqueness in the single H/V inversion procedure, using additional, complementary information is convenient. Since the dispersion curves are already computed for the theoretical H/V, the joint inversion including the fundamental Love and Rayleigh modes has no additional computational cost if data are available. As the H/V is a local measure and the dispersion curves imply at least two stations, being then non-local, the joint use of these measurements increases the constrains on the problem helping to mitigate the non-uniqueness. We show in Fig. 4 how the dispersion curves (DC) may enable in a simple case to distinguish the true solution from false ones. To account for all the information that both the H/V and the DC can supply, the cost function includes both contributions:

where \(\vartheta = \frac{n}{n + m}\) is used to equalize the new cost function. DCobs and DCth correspond, respectively, to the experimental measurement (target) and the calculated dispersion curve for the current trial model associated with the fundamental mode and evaluated at the m chosen frequencies ωi. σHV and \(\sigma_{\text{DC}}\) are the confidence intervals resulting from the computation of H/V and the dispersion curves. If they are not available, they can be replaced by a percentage of the target HVobs and/or DCobs. The idea is to minimize the cost as much as possible by exploring the universe of models using an efficient search strategy. The definition of cost function embodies the two possibilities: consider only (1) H/V or (2) the Rayleigh dispersion curves or (3) both. The non-uniqueness is always present in these types of geophysical problems. In our case, the retrieval of the underground soil velocity structure, the cost function for H/V alone has numerous minima, which are not always the same when the cost function includes only dispersion curves. Adding constrains in the cost function, in principle, reduces the velocity intervals of interest, make sharp the cost function, introducing also other minima. The challenge is to find a global minimum.

One may argue that the cost function may include Love waves or higher modes but then the observed values could be difficult to obtain. The cost function is in fact a trade-off. We use in this work the simulated annealing method but a similar optimization technique is valid as well. García-Jerez et al. (2016) and Piña-Flores et al. (2016) present a detailed and thoughtful analysis of the used inversion scheme and its implementation within a general software package.

The joint inversion was applied to the Bajo Andarax River area (SE Spain), where the sedimentary layers may produce significant effects. This area includes a part of Almería city, with about 200,000 inhabitants. Here we use data form a single station (ER20) located close to the experimental deployment. The Rayleigh DC was obtained at this site in previous studies (García-Jerez 2010; García-Jerez et al. 2010) using the SPAC technique with five concentric pentagonal setups of radius 12, 18, 25, 50, 94, and 420 m, without a central station. Ambient noise was continuously recorded for 1 h. The time windows used in correlations were 40 s long with 50% overlapping.

Figure 5 shows the result of the inversion. The fundamental Rayleigh mode associated with the profiles of Fig. 5b is very close to the experimental one. The agreement in H/V is also very good in the common frequency range. Below 0.7 Hz, the agreement is fair and some uncertainties on the depth of the profile remain. These uncertainties reflect on the gray curves associated with preliminary inversion results with an error less than two times the final error. Since three main peaks are visible on the H/V here, we assumed a profile with three layers on top on a half-space.

Fig. 5

Results of the inversion for one location at the mouth of Bajo Andarax River (SE Spain). aH/V target and theoretical for the best-fitting model. b Rayleigh (fundamental mode target and theoretical dispersion curve for the best model. c Velocity profile and gray curves are those of trial models generated by the iterative method. Misfits are less than two times the final error (after Piña-Flores et al. 2016)

We review briefly the relationship between ensemble average correlations within a diffuse field and the imaginary part of Green’s functions. We argue that auto-correlations measure directional energy densities. In fact, the imaginary part of Green’s function is proportional to the power injected into the system by a unit harmonic load, including all the reflected consequences. This is the sum of an emission episode at the source and the contribution of waves returning to the source. Various properties of the fields are the same for apparently different phenomena. Namely, the radiated energy due to the excitation of system at a point is closely related to the “projection” of the diffuse energy at such point, by means of correlations that extract the appropriate information from a ubiquitous diffuse field. One open problem is identifying and enhancing the diffuse properties of seismic ambient noise. In the meantime, the diffuse field assumption is useful.

We propose to retrieve the underground velocity structure by inverting H/V using the Simulated Annealing Method if little or no a priori information is available, or by local techniques when the misfit is close to an acceptable value, or when some information is already available. The non-uniqueness of the solutions associated with the single H/V inversion is much reduced, without additional computational cost, by doing the joint inversion of dispersion curves and H/V. The experiment for Bajo Andarax River mouth area (Almería, Spain) is a successful example of this approach.

Acknowledgements

I thank J. Piña-Flores, C. Rodríguez and N. Vera-Chávez for their multifarious help. My appreciation to U. Iturrarán-Viveros, A. García-Jerez, F. Luzón, M. Perton, T. Yokoi, and H. Wu for their constructive comments and remarks, and to E. Plata, G. Sánchez and the team of the USI—Institute of Engineering—UNAM for locating useful references. This work has partial support by AXA Research Fund and by DGAPA-UNAM under Project IN100917.

Competing interests

The author declares that he has no competing interests. He has received a salary from the Institute of Engineering of UNAM where he is Full Professor. So far, no patent has been submitted.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.