Microtremor array method using spatial autocorrelation analysis of Rayleigh-wave data

Abstract

Microtremor array measurements, and passive surface wave methods in general, have been increasingly used to non-invasively estimate shear-wave velocity structures for various purposes. The methods estimate dispersion curves and invert them for retrieving S-wave velocity profiles. This paper summarizes principles, limitations, data collection, and processing methods. It intends to enable students and practitioners to understand the principles needed to plan a microtremor array investigation, record and process the data, and evaluate the quality of investigation result. The paper focuses on the spatial autocorrelation processing method among microtremor array processing methods because of its relatively simple calculation and stable applicability.

Highlights

1. A summary of fundamental principles of calculating phase velocity from ambient noise

2. General recommendations for MAM data collection and processing using SPAC methods

3. A discussion of limitations and uncertainties in the methods

1 Introduction

During the past few decades, researchers have made considerable progress towards the development of passive surface wave methods utilizing seismic ambient noise, or microtremors. Microtremor noise as used here is background seismic noise attributable to road traffic and machinery (> 1 Hz) and coastal and meteorological sources (< 1 Hz). The wave energy propagates predominantly as surface waves, both Rayleigh waves (detectable with vertical seismic sensors) and Love waves (detectable with horizontal sensors). Although both Rayleigh- and Love-wave microtremor methods have been published, most examples (e.g., Foti et al. 2017) use vertical sensors limited to Rayleigh-wave analysis; hence, we confine discussion here to this category. We group all multi-sensor methods under the common name of microtremor array measurements (MAMs) because two-dimensional (2D) arrays are usually used for calculating phase velocity from the ambient noise. Aki (1957) first proposed the study of surface waves encoded in ambient noise and presented a theory of spatial autocorrelation, or spatially averaged coherency (SPAC).

Okada (2003) developed a large-scale MAM method based on the SPAC method in order to estimate deep shear-wave velocity (Vs) structure. Among many different methods that have been developed to retrieve surface-wave phase velocities from ambient noise, SPAC is one of the most widely used methods because of its robustness, ease of implementation, and computational efficiency. A number of review papers exist describing the theory, practical array methodology, example results of the MAM, and in particular different implementations of SPAC methods. Applications are described where the methods are used to investigate the subsurface at depth scales from a few meters to a few kilometers. Asten and Boore (2005), Boore and Asten (2008), and Stephenson et al. (2005) present comparisons by up to 14 independent groups who analyzed microtremor data in the Santa Clara Valley, California. Cornou et al. (2007) summarized results of an international blind test study of both synthetic and field data from four synthetic and four field sites. The InterPACIFIC project, funded by a consortium of European government sources, obtained detailed non-invasive passive and active seismic field measurements, plus invasive borehole studies, at three European sites representing soft sediments, firm glacial sediments, and rock (Mirandola, Italy; Grenoble, France; and Cadarache, France; respectively). In this European project, a total of 14 international groups interpreted the data, and the results were summarized in Garofalo et al. (2016a) and Garofalo et al. (2016b). Foti et al. (2017) distilled the wisdom of the multiple groups into a practical guide for future surveys. Finally, the Consortium of Strong Motion Observation Systems (COSMOS) sponsored a review paper by Asten and Hayashi (2018; hereafter AH18) and a series of blind trials where 34 international groups interpreted microtremor data from four diverse geological sites (Asten et al. 2019a; 2021a). The results of one further international blind trial were presented at the 2021 Effects of Surface Geology on Seismic Motion Symposium (Asten 2021a, b).

This document seeks to update the SPAC methodology based on experience from the COSMOS blind trials of 2018, the results of which are summarized in Section 12.

2 Principles of the MAM technique

Because the source locations generating microtremors are generally unknown, it is challenging to determine their propagation direction and thus phase velocities of surface waves composing the microtremors. Aki (1957) proposed the SPAC method in which microtremor data acquired with an array of sensors were statistically analyzed, enabling the extraction of phase velocity. The calculation of SPAC is in principle close to the general cross-correlation of ambient noise (e.g., Wapenaar 2004). SPAC can be considered the real part of cross-correlation divided by autocorrelation in the frequency domain. Further discussion and real data examples are provided in both Appendix A and AH18 (Section 10). The essential feature of the SPAC method is the calculation of the cross-correlation or cross-spectrum for the vertical-component signal from two sensors with separation r, followed by averaging of the cross-spectra around a circular array (Fig. 1). Dividing the cross-spectrum by the square root of the product of the power spectra of the two signals derives the complex coherency. For multiple pairs of stations distributed in azimuth, an average of the complex coherencies yields the spatially averaged coherency (i.e., SPAC). We note that summation around a geometric circle is not a strict requirement; we may equally use average coherencies summed for pairs of stations around the perimeter of a polygon or pairs selected from a more complex distribution of sensors such as those shown in Fig. 2, with the only provisos being similar separation distances for each station pair and a regular distribution of azimuthal angles for the pairs of stations used in the summation. Foti et al. (2017, supplement Fig. A2.2) further discuss attributes of different array geometries based on case histories.

Fig. 1

SPAC calculates cross-correlations or complex coherencies between each pair of sensors with separation r and averages these around a circle. Use of three sensors on the circle is generally regarded as a minimum. A larger number of sensors provide redundancy in case of problems with a sensor or the siting of a sensor. Where sources of the microtremor surface wave energy are highly directional, the azimuthal averaging process is improved by using an odd number of sensors on the circumference, e.g., 5, 7, or 9 sensors (Okada 2006)

Fig. 2

Examples of array geometries which have been commonly used for microtremor studies using SPAC analysis. (Modified from AH18). a Circular shape. b Nested triangular shape. c L-shape. d Sparse nested triangles (sizes increase x3). e Common-base nested triangles (sizes increase x3)

When the number of station pairs sample the wavefield sufficiently in its azimuthal distribution, the real part of SPAC (SPAC(r,ω)) becomes equal to a Bessel function (Aki 1957).

$$\mathrm{Re}\left(SPAC\left(r,\omega \right)\right)={J}_{0}\left(\frac{\omega }{c\left(\omega \right)}r\right)$$

(1)

where c(ω) is the fundamental mode phase velocity at an angular frequency ω. The left term of Eq. (1) is calculated from the vertical component of observed microtremor data, J₀ is the Bessel function (zero order, first kind), and r is the interstation distance. A more detailed overview of the SPAC method is given in both Appendix A and AH18 (Section 1).

In practice, the phase velocity is calculated by comparing the theoretical relationship in Eq. (1) with the measured SPAC and adjusting the phase velocity c(ω) to minimize the error. The velocity which minimizes this error can be considered the phase velocity at the angular frequency ω. The imaginary part of the SPAC(r,ω) ideally goes to zero in the directional average of ambient noise; the non-zero imaginary part can be used in an assessment of azimuthal distribution of the wavefield, in assessment of statistical noise in the SPAC data (Asten 2006), and as an additional control on phase velocity estimation (Cho 2020).

We can compare the SPAC Bessel functions in Eq. (1) against either the angular frequency, ω, or against the sensor separation, r (if there are multiple sensor separation distances). The SPAC against the angular frequency ω generally does not take the shape of a Bessel function because the phase velocity c(ω) is a function of frequency. Conversely, the SPAC against the sensor separation, r, always takes the shape of a Bessel function. Okada et al. (1997) named the latter plotting method extended spatial autocorrelation (ESAC or ESPAC).

SPAC, as defined in Eq. (1), can ideally only be applied to isotropic arrays (commonly where the inter-station sensor spacing is less than 10% difference among the averaged sensor pairs), such as circles or triangles, and should not be applied to anisotropic (or asymmetric) arrays, such as lines or L-shapes, due to the requirement for directional averaging. However, if we assume for a particular site that microtremors come from all directions equally over a 180° azimuthal spread, the directional average in the Eq. (1) can be calculated even if arrays are anisotropic (Capon 1973). Generally, the direction of propagation of the microtremor field is not stable. Averaging microtremor data over a sufficiently long recording time may enable us to calculate the directional average of Eq. (1) correctly even if the arrays are anisotropic, such as an L-shape or T-shape.

3 Current state of practice of MAM

SPAC has been widely applied to both engineering and scientific investigations for years (see list of references in the introduction above). The depth of investigation has varied from several meters to several kilometers in many of these published studies. The method has been used in multiple international blind test studies, and the results of the tests quantitatively demonstrate a lack of agreement, which has been attributed to both differing data acquisition and processing techniques (e.g., Asten and Boore 2005; Garofalo et al. 2016a; b; Asten et al. 2021a). An international committee (InterPACIFIC) tasked with reviewing active and passive surface wave methods summarized the data acquisition and processing of SPAC (Foti et al. 2017). This document seeks to update the SPAC methodology based on experience from the COSMOS blind trials of 2018 (results summarized in Section 12 and more fully in Asten et al. 2019a; 2021a). A wide range of different types of seismometers are commercially available, and multiple processing software packages are available, both free and for purchase (e.g., see Table 1 of Asten et al. 2021a, and detailed data descriptions in Asten et al. 2021b), which are briefly discussed in Sections 6 and 7 of this paper.

Table 1 Typical data acquisition system and investigation depth for unconsolidated or weakly consolidated sediments

Given the limitations on spatial and directional averaging of microtremor signals, as discussed in Section 2, the recommended field practice is to use isotropic arrays (also termed symmetric arrays) as a default option unless prior studies have established sufficient azimuthal spread of wave energy to provide directional averaging such that Eq. (1) is a valid assumption. Section 7.3 and Asten et al. (2019a; 2021a, Site 3) show a clear example of unsatisfactory Vs interpretation where this condition was not satisfied.

4 Limitations

Like most other geophysical investigation methods, inversion of dispersion curves obtained from SPAC is generally non-unique and includes uncertainty. Both horizontal and vertical resolution decreases proportionally as investigation depth increases, and it is in part a function of the sensor spacing, as can be seen in Eq. (1). The SPAC method assumes a one-dimensional (1D) layered structure, and the applicability to 2D structure, three-dimensional (3D) structure, or complex topography has been rarely investigated; Claprood et al. (2011; 2012) provide examples across the fault line of a rift valley in Tasmania. Chávez-García et al. (2014; 2018) studied effects on the ambient noise field of faults bounding a sedimentary basin in Greece.

The SPAC method uses ambient seismic noise, and the quality of data significantly depends on the character of the ambient noise (microtremor) wavefield at the site. The method may not work at a site with a low ambient noise level, such as a rural area or any region too remote from microtremors generated typically in urbanized areas (frequencies above ~ 1 Hz) or natural sources such as oceans (frequencies below ~ 1 Hz). Where ambient noise levels are low, a common solution is to drive a field vehicle back and forth in the vicinity of the array at locations sufficiently outside and azimuthally distributed around the array to satisfy the necessary assumptions of plane-wave and multi-directional seismic wave propagation across the array. Smith et al. (2013, their Fig. 5) provided an example using this approach. Roberts and Asten (2008) and Maranò et al. (2017) independently concluded that vehicle sources at a distance of twice the array radius from the array center provide a useable approximation to the desired plane-wave seismic wave field.

The data processing of SPAC assumes that either the receiver array or the direction of ambient noise propagation is isotropic. In practice, it is usually a combination of both array design and an azimuthal spread of the ambient noise direction of propagation giving sufficient averaging for SPAC methods to be effective. As a cautionary note, we advise that directionally anisotropic arrays (also called asymmetric arrays) such as linear arrays may not work where ambient noise propagation has strong directivity.

While isotropic 2D arrays are strongly recommended for SPAC analysis, there are situations where logistics, site access, or loss of a non-performing sensor results in an anisotropic array. Extensions to the SPAC method allowing use of arbitrarily asymmetric array data have been proposed, including the MSPAC method (Bettig et al. 2001) and the krSPAC method (Asten et al. 2019b). Section 10 gives an outline and example of the latter method.

5 Multiple modes of Rayleigh-wave propagation

Rayleigh waves can propagate over a solid or layered medium in multiple modes, each mode having a phase velocity dispersion curve that is frequency dependent. When Rayleigh waves are generated by surface sources and a layered earth has shear-wave velocity contrasts below about 2:1, the amplitude of the fundamental mode (R₀) dominates higher modes. However, where larger velocity contrasts occur, and especially when buried low-velocity layers exist beneath near-surface layers having higher Vs, first and second higher modes (R₁, R₂), or even higher modes, may carry a dominant proportion of wave energy for some frequency bands (e.g., Claprood et al. 2011; Ikeda et al. 2012; Salloum et al. 2014; Hayashi 2019).

The partition of energy between modes is dependent on three factors: the layer velocities (i.e., specific acoustic impedance contrasts), the attenuation, and the nature and depth of the wave sources (whether vertical impacts or other types of sources). If we assume that sources are from vertical impacts at the surface of the Earth and that the Earth consists of laterally uniform elastic homogeneous layers, the energy partition between different modes can be theoretically computed (e.g., Ikeda et al. 2012; Hayashi and Craig 2017), and we can define an “effective mode,” often labeled R_e. Figure 3 shows the theoretical first three Rayleigh-mode phase velocity dispersion curves (solid lines), their relative amplitude (dashed lines), and the dispersion curve for the effective mode for a high-contrast interface at 20-m depth. For most of the frequency bands of interest, the R₀ mode dominates, but between 5 and 10 Hz a significant fraction of the Rayleigh-wave energy propagates in the R₁ mode, with the result that the R_e dispersion curve separates from R₀ and shows an upward notch to higher velocities in that frequency band. It implies that an interpretation taking only the fundamental mode into account may provide incorrect velocity-depth profiles. In many field studies over high-contrast layer boundaries, observed dispersion curves show the same notched shape associated with frequency-dependent changes in energy partition between modes.

Fig. 3

a Example of theoretical dispersion curves for three Rayleigh modes over a high-contrast interface. Solid and dashed red, yellow, and green lines show dispersion curves and their relative amplitudes for modes R₀, R₁, and R₂, respectively. Blue line is effective mode R_e (e.g., Ikeda et al. 2012; Hayashi and Craig 2017). b Model has two layers, Vs = 400 m/s, thickness 20 m, overlying a hard basement Vs = 2000 m/s. From Foti et al. (2017, Appendix 5)

The consequence of such energy (power) partitioning between modes is typically small for studies of geological sections with small shear-wave velocity contrasts across layer boundaries (i.e., less than about 2:1). Interpretation in such cases may be achieved considering only the fundamental mode of propagation. However, when larger shear-wave velocity layer contrasts exist (> 2:1), correct interpretation may require analysis in terms of the “effective mode” dispersion curve (blue solid line in Fig. 3) in which the higher modes are explicitly modeled. Such analysis is not perfect since it assumes the power partition between modes follows a theoretical model for vertically acting surface sources on a laterally homogeneous earth. However, practical experience indicates it is a useful approximation when strong layer contrasts are suspected. Other exceptions may occur where energy sources are not at the surface (e.g., buried machinery) or when lateral changes in Vs produce mode conversions in Rayleigh waves. Foti et al. (2017, Appendix 5) provide an extensive discussion of these issues.

6 Data collection

6.1 Equipment

Two types of data acquisition systems are mainly used for obtaining microtremor array data: (1) a multi-channel seismograph system, where dozens of receivers typically are connected through multi-channel cables (wires); and (2) a wireless array of self-contained seismographs that includes a GPS clock with each sensor so that all units can be synchronized over any distance. The multi-channel seismographs generally use receivers which are velocity meters (geophones) with a natural frequency of 1 to 4.5 Hz. Wireless seismograph systems commonly use either accelerometers or velocity sensors. The multi-channel seismographs are mainly used for smaller (< 100 m) arrays, and the wireless seismographs are mainly used for larger arrays (> 100 m). Since the wireless seismographs are generally used for larger arrays that require low-frequency phase velocities, broadband sensors or geophones with a low resonance frequency are used as receivers. Multi-channel seismographs have been used for exploration seismology, and many different types are commercially available. Wireless seismograph systems are becoming popular since instrumentation has become smaller and less expensive, owing to the progress of electronics. These newer systems are commonly referred to as nodal arrays. The growing availability of wireless seismographs with broadband sensors enables us to perform SPAC measurements efficiently with large seismic arrays. The choice of acquisition systems strongly depends on the noise level, geology of sites, and investigation depth. Table 1 summarizes typical data acquisition systems and achievable investigation depths over unconsolidated or weakly consolidated sediments. It is important to choose appropriate equipment and verify suitable data quality during fieldwork immediately after data collection.

6.2 Deployment

Our recommendation for passive surface wave methods is to use 2D arrays such as those shown in Fig. 2 or more generally by Foti et al. (2017) and AH18 because SPAC essentially depends on directional averaging of the wave field, as mentioned in Section 2. The circular array gives the most effective azimuthal averaging, especially where an odd number of stations are used on the circle circumference (Okada 2006). Poggi et al. (2017) give examples of arrays consisting of multiple circular arrays applied to challenging stiff soil sites in the Alps. For more general applications, the sparse nested or common-base triangles are the most efficient geometry since these give sufficient azimuthal averaging in most cases. Zhang and Pankow (2021) give examples of how square sub-arrays extracted from a larger grid of stations allow the SPAC method to be applied over a region, with the result that a quasi-3D geological model can be constructed from the set of 1D interpretations from each sub-array. In areas of restricted access, L-shaped arrays can also be useful. In some cases, useful SPAC data can be obtained with a linear array or even a two-station array (Hayashi 2009). However, it should be noted that linear arrays require the propagation of ambient noise to be omni-directional over 180° of azimuth, and that 2D regular arrays are generally desirable.

Array size (or maximum receiver spacing) and the maximum wavelength obtained from SPAC are directly related. As a rule of thumb, the maximum wavelength is generally between 2 and 4 times the array size regardless of array size and shape (Cornou et al. 2007; Hayashi 2019). It is generally agreed that the maximum penetration depth of surface wave methods is roughly in the range of 1/2 to 1/4 of maximum wavelength (e.g., Xia et al. 1999). We therefore conclude that the approximate penetration depth of the passive surface wave methods using SPAC is between 1 and 1.5 of the maximum receiver separation for a given array (Foti et al. 2017). It should be noted that the penetration depth of SPAC strongly depends on the noise level and geology of sites. The penetration depth might be much shallower compared with the maximum receiver spacing at low noise level sites. The use of broadband sensors and long record length are necessary to assure low-frequency data and a large penetration depth.

The minimum penetration depth (or near-surface resolution) achievable with SPAC is dependent on the shortest wavelength resolvable. The short wavelength limit is in part dependent on the useable spectral bandwidth in the Rayleigh-wave data, which in turn is dependent on the choice of the SPAC software processing algorithm used (Asten et al. 2021a).

We propose the guidelines for minimum and maximum depths of investigation D_min, D_max, to be

$$\begin{array}{c}{D}_{\mathrm{min}}=\lambda min/3=V/\left(3\cdot fmax\right),\\ {D}_{\mathrm{max}}=\neg \lambda max/2=V/\left(2\cdot fmin\right)\end{array}$$

(2)

where λmin and λmax are the minimum and maximum wavelengths, fmin and fmax are the minimum and maximum frequencies provided by SPAC analysis respectively, and V is the phase velocity observed at a given fmin or fmax. Further discussion is provided in AH18 Section 3.

6.3 Recording

Longer recording time (and thus more data) is better for statistical analysis in SPAC. However, the need for long data acquisition times decreases the field efficiency of the method; therefore, recording microtremor data with an appropriate data length is a very important factor in practical survey design. Several experimental studies (i.e., Hayashi et al. 2006; Hayashi 2009; 2019) demonstrate that 10 to 20 min of data are usually enough for arrays with a spacing of less than 100 m. The record length needs to be increased for larger array sizes, although data quality varies between sites. Typical recording times are 30 to 40 min for an array size of several hundred meters, and a couple of hours for arrays larger than 1 km. Data longer than several hours usually does not improve data quality significantly; provided sources of microtremor energy are constant.

However, some surveys such as at urban sites can have significant problems with local seismic sources (e.g., machinery) existing within an array; these can lead to near-field nonplanar wave propagation resulting in low interstation coherencies. Figure 4 shows an example where nighttime recording is strongly preferable due to absence of local daytime sources. Another instance where long recordings may be desirable is at low-noise remote sites where energy sources for microtremor energy are sparse. Schramm et al. (2012) give such an example where the most useful microtremor energy in the desert environment over 24 h was from late afternoon recordings; that result may have been associated with increased wind energy late in the day. These examples illustrate exceptional situations where recording overnight may be advantageous if site logistics can be arranged.

Fig. 4

Comparison of observed coherencies obtained during daytime and nighttime in an urban area (downtown Seattle, WA, USA), using receiver distances of a 150 m and b 200 m. Coherencies observed in the nighttime (blue lines) for low frequencies (< 1.5 Hz) are clearly larger than those observed in the daytime (orange lines)

7 Data processing

7.1 Phase velocity calculation from microtremor waveforms

Two approaches are accepted for data processing steps in SPAC analysis. In the first approach, the recorded vertical-component microtremor time-series data are divided into multiple time blocks (typically a couple of minutes) that can be overlapped by a small amount of time, or by a proportion of the time block length. The time series of all stations must be aligned to a common time origin before selecting the time blocks. As a quality-control step, selection of these blocks will typically exclude bursts of non-stationary noise, such as those from local moving vehicles. A fast Fourier transform (FFT) is applied to each time block to transform the time domain waveform data to the frequency domain. Coherencies between pairs of sensors are first calculated for each time block, and then coherencies of all time blocks are averaged for a set of sensor spacings (Appendix A). If the coherencies are averaged over many time blocks or over a sufficiently long time and over many different azimuthal pairs, it can be considered to approximate a Bessel function, as shown by Eq. (1). The processing can be done in the time domain, although most of the processing algorithms currently available use coherencies computed in the frequency domain. A more detailed overview of the SPAC method is given in Appendix A.

In the second approach, after selection of a single long block of data (with deleterious features such as local impulsive spikes removed), perform the FFT on the single block; compute the coherency spectra for pairs of stations; and average the complex coherencies over the specified frequency windows (Asten 2006). The two approaches are equivalent.

Figure 5a shows an example of coherencies as a function of interstation sensor distance. A phase velocity can be determined at each frequency so that the difference between both sides of Eq. (1) is minimized. Figure 5b is a phase velocity image in the frequency domain that shows the error between observed coherencies and theoretical Bessel functions. This series of phase velocities defines a dispersion curve (phase velocity-frequency observations; Fig. 5c). A velocity model for shear-wave velocity (Vs) versus depth is obtained by inversion, as described in the next sub-section.

Fig. 5

Outline of processing based on SPAC. a Coherencies against distance (separation between receivers). b Phase velocity image in frequency domain. Differences of color indicate differences of error between observed coherences and theoretical Bessel functions shown in equation Eq. (1). Blue color indicates small error and red color indicates large error. Red dots show phase velocities that give minimum difference between both sides of Eq. (1). c Dispersion curve. d Initial velocity model together with 1/3 of wavelength (red circles) converted from observed phase velocities

7.2 Inversion

A Vs model is generally obtained from dispersion curves (first calculated through the waveform processing outlined in the previous section), by a non-linear inversion. Most inversions require some initial models or search areas to estimate appropriate models. These models and search areas are generally created by a simple wavelength transformation (Xia et al. 1999), in which wavelengths are calculated from phase velocity and frequency divided by 3 and plotted as depth. Figure 5d shows an example of an initial velocity model overlain by the 1/3 wavelength-phase velocity pairs calculated from the dispersion data in Fig. 5c.

There are many non-linear inversion methods, with iterative non-linear least squares, genetic algorithms, and simulated annealing among the most prominent (e.g., Menke 1984). All inversion methods try to reduce the error between the observed and predicted phase velocities.

$${\sum }_{i}^{N}{\left({f}_{i}^{\mathrm{obs}}-{f}_{i}^{\mathrm{cal}}\left(V{s}_{1},V{s}_{2},\cdot \cdot \cdot \cdot V{s}_{M}\right)\right)}^{2}={\sum }_{i}^{N}{\left({f}_{i}^{\mathrm{obs}}-{f}_{i}^{\mathrm{cal}}\left(x\right)\right)}^{2}\to \mathrm{Minimum}$$

(3)

where Vs₁, Vs₂, ・・・Vs_M are the Vs values for the 1st, 2nd, and Mth layers, respectively. N is the number of observed phase velocity data; f^obs represents phase velocities obtained from observed waveform data; and f^cal is the set of fundamental-mode theoretical phase velocities for the Vs model. In the non-linear least squares method, the iterative process changes the Vs until a good fit is obtained between the observed and calculated phase velocities. Theoretical dispersion curves are generally calculated by a matrix method (e.g., Saito and Kabaswa 1993). Since this is the most commonly used method, we summarize the details of the iterative non-linear least squares method in Appendix B.

An alternative approach to fitting observed and dispersion curves is to fit observed and predicted SPAC curves directly. As with the previous method, the model SPAC curve requires computation of a model dispersion curve for a specified layered earth, but it has some advantages in avoiding the highly non-linear step of obtaining an observed phase velocity dispersion curve from a SPAC spectrum. Known as multi-mode SPAC, MMSPAC, or direct-fitting SPAC, the approach may have an advantage in inverting Eq. (1) when, by definition, the inverse Bessel function is multi-valued (Asten et al. 2004; Asten 2006). Figure 6 illustrates the alternative processing and inversion paths used by ESAC or SPAC dispersion curve and MMSPAC direct-fitting methods. Zhang et al. (2019) provides a detailed discussion of inversion methodology for the direct fitting method. AH18 Section 5 provides a more detailed discussion of both processing paths. Section 9 provides field examples demonstrating both approaches.

Fig. 6

Alternative processing and inversion paths used by ESAC or SPAC dispersion curve methods (at left) and by MMSPAC direct-fitting methods (at right). A triangular array (top right) provides SPAC spectra for two station separations r1 and r2, which can be inverted or fitted simultaneously using either of the two processing paths. (Modified from AH18)

The calculation of theoretical dispersion curves requires P-wave velocity (Vp) and density in addition to Vs. Since the effects of Vp and density in the dispersion calculation are much smaller compared to Vs (Xia et al. 1999), Vp and density are commonly automatically changed with Vs or set to reasonable constant values. In the absence of other information, Vp above groundwater level is often assumed to be double Vs, and beneath the groundwater level Vp may be set to reasonable constant values, or it may be automatically changed based on an empirical relationship (e.g., Kitsunezaki et al. 1990; AH18 Section 3). Vp-Vs relationship shown in Kitsunezaki et al. (1990) can be written as,

$$Vp=1.11\cdot Vs+1290\left(\mathrm{m}/\mathrm{s}\right)$$

(4)

Similarly, density may be assumed, or it may be automatically changed based on empirical relationships (e.g., Ludwig et al. 1970).

The choice of method for Vp and density is not critical because dispersion curves are generally insensitive to these parameters. An exception applies at the water table; because Vp above and below the groundwater level shows a large change in unconsolidated soil, it follows that the location of the water table may have a large effect on an estimated Vs model in soft sediments. In such cases, it is desirable to obtain reliable groundwater information from other investigations, such as borings. If there is no reliable groundwater information, it may be better to assume the water level to be at the ground surface. Shallow groundwater implies high shallow Vp, which results in interpretations of lower Vs; these lower values can be considered conservative estimates from an engineering point of view.

Most active and passive surface wave methods assume that dispersion curves are dominated by the fundamental mode (R₀). The higher modes may, however, dominate in several types of velocity structures, such as in situations where a high-velocity layer overlies a low-velocity layer, or where high-velocity layers are embedded between low-velocity layers. It is therefore necessary to take the higher modes into account in the inversion of dispersion curves in these cases (see previous section).

Analysis including higher modes remains an active research topic. In SPAC analysis, the effective mode (R_e) introduced in Section 5 is often used to account for higher modes (e.g., Obuchi et al. 2004; Ikeda et al. 2012; Foti et al. 2017 Appendix 5). The effective mode is the weighted average of the fundamental and higher modes. SPAC analysis including higher modes can be written using the effective mode of surface waves C_s(ω).

$${J}_{0}\left\{\frac{\omega r}{{C}_{s}\left(\omega \right)}\right\}=\frac{{\sum }_{i=0}^{n}{p}_{i}\left(\omega \right)}{P\left(\omega \right)}{J}_{0}\left\{\frac{\omega r}{{C}_{i}\left(\omega \right)}\right\}$$

(5)

where p_i(ω) is the power fraction of the ith mode of the surface waves, n is the number of modes, and ω is the angular frequency. P(ω) is the sum of power of all modes. C_i(ω) is the phase velocity of the ith mode of surface waves, and r is the receiver spacing.

To calculate the effective mode phase velocity (C_s(ω)), medium response and phase velocity (C_i(ω)) are calculated by the matrix method. The medium response can be considered relative amplitude among fundamental and higher modes of surface waves used as the power (p_i(ω)) in Eq. (5). In SPAC processing, observed coherencies, or the left-hand side of Eq. (1), are calculated from ambient noise, and the effective mode phase velocity (C_s(ω)) is used to calculate the Bessel function (the right-hand side of Eq. 1) instead of the fundamental mode phase velocity (c(ω)). Figure 7 shows an example of a dispersion curve including higher modes and an estimated velocity model. The fundamental mode and higher modes up to the 5th mode were considered in the calculation of the theoretical dispersion curves. Observed phase velocities are shown as white circles with a solid red line, and the modeled effective mode R_e phase velocities are shown as yellow circles. Solid and broken lines indicate fundamental and higher mode theoretical dispersion curves and their relative amplitude (response of the medium). We can see that a bump on the dispersion curve around a frequency of 1 Hz is due to the relatively large amplitude of the first higher mode. We cannot simply apply Eq. (5) to the SPAC data from different receiver spacings. Instead, a direct-fitting method described below enables us to compare coherencies and Bessel functions for different receiver spacings individually. Additional examples of fitting observed and model SPAC curves using assumed R₀ and R_e propagation are provided in Sections 9 and 10, and in AH18 (their Figs. 4 and 7).

Fig. 7

a Example of an observed dispersion curve (solid red line) together with model (effective mode R_e) dispersion curves (yellow circles) computed using the Rayleigh-wave fundamental and five higher modes. Solid lines indicate fundamental and higher mode model dispersion curves, and broken lines indicate relative medium responses (e.g., Ikeda et al. 2012; Hayashi and Craig 2017). b The interpreted shear-wave velocity model together with the 1/3 wavelength calculations (red circles) converted from observed phase velocities

The iterative non-linear least squares method is computationally fast and generally stable if appropriate initial models were provided. It should be noted that the method can be very easily trapped in a local minimum. The method also requires the calculation of partial differentials. Fortunately, the fundamental mode dispersion is generally continuous, and the partial differentials can be calculated. The effective mode or dispersion curves including the higher modes, however, are often discontinuous, and the calculated partial differentials may not be meaningful. In such a case, the iterative non-linear least square method cannot be applied. Global search methods, such as a Monte Carlo method (e.g., Saifuddin et al. 2018), a genetic algorithm (e.g., Yamanaka and Ishida 1996), or a simulated annealing algorithm, are often used instead of the iterative non-linear least squares method to avoid being trapped in a local minimum. The methods are computationally slow but do not require the calculation of partial differentials and can be applied to dispersion curves including higher modes (Hayashi 2019). The various inversion methods all have strengths and weaknesses, and we recommend applying several different inversion methods with different initial models, search areas, constraints, and parameters. Among various inversion methods, we prefer the non-linear least squares method with fundamental mode of dispersion curve as the first choice since the method is relatively simple, fast, and stable. The genetic algorithm or simulated annealing with a higher mode dispersion curve may be the second choice if field data are found to include significant higher modes. MMSPAC or direct fitting SPAC can be also an alternative for taking the higher modes into account.

7.3 Uncertainty analysis

The inversion of dispersion curves obtained from SPAC is generally non-unique, like most other geophysical methods. In this section, we mainly summarize the causes and consequences of uncertainty associated with the calculation of coherencies or dispersion curves from ambient noise by SPAC methods. As mentioned earlier, SPAC requires either omni-directional noise propagation or receiver array geometry to be isotropic. Asymmetric arrays, such as a linear array, may not work where noise propagation has strong directivity. The 2D arrays shown in Fig. 5 generally provide reliable dispersion curves even if the noise propagation has strong directivity. Figure 8 shows an example of coherencies obtained from ambient noise with strong directivity. An equilateral triangle with four receivers is used in the data acquisition. Three coherencies are calculated between each corner and center receiver. In the ideal case, the three coherencies should be identical if the noise is isotropic. However, we see that the three coherencies have large differences, which indicates strong noise directivity. To avoid problems caused by such directivity, the use of isotropic arrays, such as circular or triangular shape shown in Fig. 2, is strongly recommended. If site access is such that only a linear array can be used, measuring with several different directions, with different time periods, or by comparing with active methods, may increase the reliability of the investigation.

Fig. 8

Example of coherencies obtained from the noise with strong directivity (site 3 of Asten et al. 2021a). a Array used in data acquisition. b Coherency spectra calculated between individual pairs of receivers. It appears that the ambient noise mainly propagates from bottom left to top right

Inappropriate array size or receiver spacing may decrease accuracy and reliability. Insufficient array size particularly tends to result in phase velocity estimates that are systematically underestimated at low frequencies or long wavelengths. Figure 9 shows the example of dispersion curves obtained from three different array sizes at the same site. We can see that the phase velocity tends to be biased to low values when the wavelength is longer than approximately double the array size. To avoid the problem, the use of a sufficiently large array is crucial. It is strongly recommended to use a maximum receiver spacing larger than the depth of investigation. Sufficient ambient noise level in the frequency range of interest and the use of appropriate sensors are also important.

Fig. 9

Example of dispersion curves (a) obtained from three different size of arrays (b) at the same site. Solid lines in a indicate wavelength approximately double array sizes. Note that apparent phase velocities where the wavelength is longer than 1.5 to 2 times the array size are biased by noise to low values; these must be discarded before inversion to a velocity model

SPAC assumes all microtremor sources are far enough away from an array of measurements so that plane wave propagation can be assumed. Sources placed within or near the array may cause serious problems in phase velocity calculations. Traffic noise in urban areas is a typical example of such undesirable sources. The traffic noise is generally large during daytime and small during nighttime; therefore, it may be better to perform data acquisition of larger arrays at nighttime in urban areas. Figure 4 shows a comparison of observed coherencies obtained during daytime and nighttime in an urban area, and in this example, we can see that coherencies during the daytime are reduced to unusable small values for frequencies below 1.5 Hz. By contrast, coherencies from data acquired during nighttime are useable down to a frequency of 0.2 Hz.

8 Alternative data collection and processing methods

There are many different data acquisition and processing methods to estimate Vs structures from seismic ambient noise (Table 2). Socco and Strobbia (2004) describe various modalities for the processing of seismic data in terms of surface wave analysis (like tau-p, frequency decomposition, omega-c, and frequency wavenumber). People have used SPAC and frequency wavenumber beamforming (FK) (Capon 1969) for many years. The most important advantage of FK processing over other methods is that the method can estimate propagation direction of ambient noise and can separate multiple modes of Rayleigh-wave propagation without requiring the concept of the effective mode. An important disadvantage of FK processing is that the method requires a relatively large number of sensors (at least 7, 10 to 14 preferred) compared with SPAC. Poggi et al. (2017) gave examples of the detailed study of the wave-field possible when using arrays with 19 to 30 sensors. Wathelet et al. (2018) also studied the three-component FK method with both synthetic and field data and showed its advantages in identifying multiple Rayleigh modes. AH18 (Section 9) discussed the strengths and weaknesses of the SPAC and FK methods. Seismic interferometry (SI), proposed by Wapenaar (2004), has gained popularity to calculate phase or group velocity from ambient noise array data. The SI approach allows users to calculate group velocity from data recorded by spatially aliased receiver arrays, while other methods generally require a spatially unaliased receiver array. A clear disadvantage of the SI method is that it requires a much longer record length of microtremors compared with SPAC. Appendix A and AH18 (Section 10) discussed similarities and differences between SPAC and SI. Cho et al. (2013) proposed a centerless circular array (CCA) method and claimed that the method could estimate phase velocities for wavelengths up to 100 times the diameter of the circular array. Louie (2001) proposed a passive surface wave method using a linear array, named the refraction microtremor (ReMi) method (see paper by Louie et al. 2021, this volume). This method is widely used to estimate average Vs down to 30-m depth (Vs₃₀) for engineering purposes since the method is relatively simple and easy to carry out using conventional multi-channel seismographs.

Table 2 Data collection and processing methods to estimate Vs models from seismic ambient noise

The use of Love-wave dispersion curves, obtained from the horizontal component of ambient noise, has been proposed as a method to decrease the uncertainty of Vs inversion and increase the accuracy of a given investigation. Although numerous attempts have been made by researchers to extract Love wave dispersion curves from three-component ambient noise (e.g., Aki 1957; Yamamoto 2000; Okada 2003), the use of horizontal-component data still remains in the research and development phase. The recorded vertical component of ambient noise only consists of Rayleigh waves, whereas horizontal components consist of both Rayleigh and Love waves. It is thus necessary to estimate the Rayleigh and Love amplitude ratio to calculate Love-wave velocity, and this requirement makes the use of horizontal-component data in dispersion-curve analysis difficult.

All the data correction and processing methods summarized above have strengths and weaknesses, and blind tests such as that summarized in Section 11 suggest no microtremor array data processing method is wholly superior over other methods. The best method depends on site conditions, investigation purposes, and available resources, among other factors. We recommend applying several different processing methods to the same data and comparing the resulting dispersion curves. In this way, the uncertainty in calculating the dispersion curve can be decreased, and the reliability of the resultant Vs structures can be increased. It is important to choose array configurations to which several processing methods can be applied. Among various processing methods, we prefer SPAC as the first choice since the method is generally robust and stable and works with a small number of sensors. FK beamforming may be the second choice if we need to estimate the direction of ambient noise propagation and if many sensors are available. SI is a useful alternative if the use of wide sensor spacings has resulted in spatially aliased data.

9 Case studies

A large number of case histories exist where MAM investigations are supplemented with three-component horizontal-to-vertical spectral ratio data (HVSR). The HVSR method is reviewed in detail by Molnar et al. (2018; 2021). We include HVSR data for the Mexico City case history discussed in this section.

The M_W 8.0 earthquake that struck Mexico City on 19 September 1985 caused severe damage although the city was 400 km from the epicenter. The main reason for this is that the city is located on a very soft sedimentary basin (Abe 1986). The distribution of these soft sediments has been delineated by drillhole records and microtremor measurements (Lermo and Chávez-García 1994). Vs profiles were measured by downhole seismic logging (Seo 1986) at only seven sites, and the regional Vs structure of the whole basin was not delineated. To delineate the Vs structure of the basin below central Mexico City down to a depth of approximately 200 m, three-component microtremor data were analyzed using both SPAC and HVSR (Hayashi et al. 2011) at six sites. The MAM surveys were conducted with 25- to 650-m equilateral triangular arrays (Fig. 10a). Figure 10b to d show an example of a dispersion curve, an HVSR spectrum, and a Vs model obtained in the center of the Mexico City basin. A clear dispersion curve was obtained between 0.3 and 7 Hz (Fig. 10b), and a HVSR peak was observed at a period of 3.5 s (Fig. 10c). A joint inversion of the dispersion curve and HVSR spectrum was applied to observed data, and a Vs model down to a depth of 200 m was obtained (Fig. 10d). Low-velocity layers with Vs approximately 30 and 100 m/s existed to depths of approximately 20 and 80 m, respectively. It seems that the peak period of 3.5 s in the HVSR spectra is due to the two shallow low-velocity layers with Vs of 30 and 100 m/s, respectively.

Fig. 10

Shear-wave velocity investigation in the Mexico City basin, using a passive seismic array. a An example of the array, b observed and theoretical dispersion curves, c observed and theoretical HVSR, and d an estimated Vs velocity model. A clear peak period of 3.5 s in the HVSR spectra is due to the two shallow low-velocity layers with Vs 30 m/s and 100 m/s. (Modified from Hayashi et al. 2011)

The direct fitting method of SPAC interpretation is illustrated using samples of microtremor data from Hollister Municipal Airport, California. The site and early microtremor surveys along with borehole data were described by Liu et al. (2000). This paper uses data acquired in a later survey with pairs of nested triangles (see example in Fig. 2b), where the triangle side lengths range from 12.5 to 100 m, supplemented by asymmetric linear arrays placed along runway margins. Figure 11a shows the location of triangular arrays at Hollister Airport, and Fig. 11 shows as an example the observed and model coherency spectra acquired with the triangle of side length 50 m; it is evident that useful data exist over a bandwidth of 2 to 17 Hz for spatial averaging of coherencies over the three radial separations (r = 28.8 m), and 2 to 12 Hz for averaging over the side lengths (r = 50 m). Section 11 illustrates how this methodology may be extended to cases where the arrays are not symmetric.

Fig. 11

a Locality for MAM triangular arrays at Hollister Municipal Airport, California. Center of the array is at 36.888543°N 121.413778°W. Yellow pins show four nested triangles. Blue shading shows one triangle side length 50 m used for this figure. b Vs profile interpreted from all arrays, including larger arrays shown in Fig. 12, using geological data provided in Liu et al (2000). c, d SPAC spectra for triangle radii and sides with interstation distances r = 28.8 and 50 m, respectively. Black line—observed SPAC; red and blue lines—model SPAC spectra for the fundamental Rayleigh mode R₀ and the effective mode R_e. The fitting is performed by least squares using the observed and model R_e curves. Thick black line—the frequency range used in the curve fitting. The radial distance recovers useful data over a frequency band 2 to 17 Hz. Side length distance likewise over band 2 to 12 Hz

Further case studies integrating use of SPAC methods and HVSR spectra are provided in AH18 (Sections 7 to 9) and elsewhere in the literature (e.g., Stephenson et al. 2009; Suzuki and Yamanaka 2010; Zor et al. 2010; Hayashi et al. 2011; Asten et al. 2014; Hayashi and Craig 2017; Chieppa et al. 2020).

The challenge of using passive microtremor methods to correctly detect and quantify low-velocity layers (LVL) situated below overlying higher-velocity sediments is difficult. Surface-wave models are subject to non-uniqueness limitations; the blind trial reported by Garofalo et al. (2016a) included one site (Grenoble, glacial sediments) having a 10-m-thick layer of soft clay under 25 m of surficial sands and gravels. The blind trial showed only four of 14 groups correctly identified presence of the LVL. Roberts and Asten (2004) describe a site where an LVL of 25 m of Quaternary alluvial deposits beneath 12 m of Quaternary basalts was resolved by SPAC methods. Farrugia et al. (2016) quantified thicknesses of 20 to 50 m of soft clays below 10 to 50 m of hard surficial limestones. Teague et al. (2018) resolved multiple layers of sand interbedded with gravels over a total thickness of 150 m, aided by geological control. AH18 contains additional discussion of the LVL challenge for microtremor methods.

10 Advanced approaches

Where logistics and instrument availability permit, dense arrays comprising circles of stations with five or more stations per circle are highly effective. Array processing using FK methods gives precise results with such arrays and has the advantage that in the presence of directional seismic noise and/or departures from symmetry in the array geometry the MAM method remains robust. Foti et al. (2017) provide a tutorial explanation of FK methods, and a set of detailed examples is given in Poggi et al. (2017).

With less dense arrays, typically used where logistical and instrumental constraints place limits on the complexity of array geometry, both the above advantages cease to exist. In addition, site access or corruption of the microtremor signal by local noise in the vicinity of a sensor may reduce a nested triangle array to an asymmetrical shape. In order to gain useful SPAC spectra from such arrays, a variant termed krSPAC (Asten et al. 2019b; Stephenson et al. 2019; Asten 2021a, b) may be used. The concept behind the krSPAC approach is that if the frequency axis of coherency curves is transformed to the dimensionless scale kr (where k is wavenumber and r is separation of pairs of stations), coherency curves can be averaged for multiple pairs of stations regardless of whether the pairs have different spacings. Details of the method are given in Asten et al. (2019b) and Stephenson et al. (2019).

Figure 12 shows how the krSPAC approach allows direct fitting to be used on an asymmetric array produced by the logistical restriction of access to runway margins at the Hollister Airport (Fig. 12a). Acquisition instrumentation was a set of 2-Hz geophones wired to a multichannel seismograph. Apart from geometrical constraints, data at many stations are partially corrupted by noise presumed to be from industrial machinery or signaling, dominant at frequencies 1.3 and 1.7 Hz. After review of noise characteristics at various stations, we selected an asymmetric four-station triangle as shown in Fig. 12a, having unequal “radii” (292, 156, 151 m) and unequal side lengths (364, 307, 295 m) to develop the Vs model shown in Fig. 12b. Figure 12c and e show the conventional SPAC direct fitting, with the layered earth model being very close to one of the models hypothesized in Fig. 11 (curve 5) in Liu et al. (2000). The spatial averaging used in the conventional method fails due to the unequal station spacing, for frequencies above about 2 Hz, giving meaningless results at these higher frequencies. However, when spatial averaging is performed in the kr domain, meaningful fitting of observed data to model data is achieved over frequency ranges up to 4 and 3 Hz, respectively, for spatial averaging over the three notional “radii” and over the three side lengths (Fig. 12d, f).

Fig. 12

a Locality for large asymmetric array at Hollister Municipal Airport, California. Yellow pins show geophone positions. Blue asymmetric “triangle” shows a four-station triangle selected for krSPAC interpretation. Yellow arrow shows the station selected as a notional “center” for purposes of constructing unequal “radial” distances to each vertex of the blue triangle. b Vs profile interpreted using data from Liu et al (2000) and this study (see also Fig. 11b). c, e Conventional SPAC spectra for asymmetric blue triangle; averaged “radial” and averaged side length distances r_ave = 200, 322 m respectively. The large variation in inter-station distances used in the spatial averages makes the conventional SPAC plots of dubious value for frequencies above 2 Hz. Black line—observed SPAC; red and blue lines—model SPAC spectra for the fundamental Rayleigh mode R₀ and the effective mode R_e. The fitting is performed by least squares using the observed and model R_e curves. d, f SPAC plots obtained using the krSPAC algorithm where coherencies for pairs of unequally spaced stations are averaged in dimensionless wavenumber space kr rather than in frequency space (as described in Asten et al., 2018). Annotations above the kr axis are equivalent frequencies. Thick black line—the frequency range used in the curve fitting of observed SPAC and modeled SPAC computed from R_e. The “radial” distances d recover useful data over a frequency band 0.3 to 4 Hz. The side length distances f give useful data over the band 0.3 to 3 Hz. The standard deviation of observed and model R_e data improves from 0.14 for conventional SPAC in e to 0.10 for krSPAC in f

11 Lessons from the COSMOS blind trials 2018–2020

The COSMOS blind trials of MAM methods involved 34 groups of analysts using about 13 different software packages to perform blind interpretations on MAM array data from four different sites in the USA, Canada, and Europe (Asten et al. 2019a; 2021a). For each site, the array data were released to analysts in four phases, allowing objective analysis of the usefulness of simple two-station arrays, sparse arrays (simple triangles), and dense arrays (circles or nested triangles). The final phase was reinterpretation and reconciliation of the geophysical data with all available geological information. The phase 1 results illustrated the dangers of erroneous interpretations of Vs when using two-station arrays in the presence of directional seismic noise. Our conclusion is that we recommend in Section 6.2 against use of such two-station arrays if it can be avoided.

The phase 2 trials (sparse triangular arrays) showed that a majority of analysts were able to get acceptable Vs profiles where the goal was time-averaged Vs data (V_S10, V_S30, or V_S100) for the site as required for general earthquake hazard codes (e.g., BSSC, 2003). However, such data were not sufficient where depth to major interfaces (e.g., base of surficial soft sediments or depth to basement rock) are required.

The phase 3 trials using dense 2D arrays demonstrated that accurate interpretation of depth to interfaces from MAM surveys requires the use of dense or multiple arrays and the use of three-component data to allow HVSR data to be used in conjunction with phase velocity dispersion data. Section 9 and Fig. 9 of this paper illustrate this process.

The blind trials did not identify any MAM software package as clearly superior; for example, the top six analysts assessing time-averaged Vs data used five different software systems. Analysts’ practice and skill were deemed important. To facilitate training in these skills, the blind trial data are available in the four phases and thus are a tool for future education in MAM interpretation (Asten et al. 2021b).

Opinion among analysts was split based on the question of whether active surface wave methods should be employed in order to define the Vs profile for near surface (e.g., upper ~ 10 m) of topsoil. Near surface definition requires high-frequency Rayleigh-wave dispersion data, typically above 10 Hz. The blind trials demonstrated that the widest useable bandwidths of microtremor data were obtained by analysts using the ESAC method and with the MMSPAC method based on direct fitting of coherency spectra; it is arguable that maximizing the useable bandwidth in passive seismic methods up to 20 ~ 30 Hz removes the need for supplementary active surface-wave measurements (e.g., Hayashi et al. 2016).

The last of these points is illustrated in the present paper. Generally, with low-noise data and an absence of lateral variations in geology, the observed SPAC curve can be interpreted to a high-frequency limit of about the third minimum of the plot (see e.g., Fig. 11c). This minimum corresponds to a wavelength of order 0.4r, where r is the minimum sensor separation in the array.

12 Conclusions

This paper summarized fundamental theory of the microtremor array method and provided general guidelines for data acquisition and processing. The results of several blind tests have demonstrated that MAM can provide reliable investigation results applicable for engineering purposes. We recommend that surveys should ideally use isotropic arrays (e.g., circular or triangular), or at least 2D arrays (e.g., L-shape), with the maximum receiver spacing being larger than the desired investigation depth. Various phase velocity calculation and inversion methods have strengths and weaknesses, and no method is conclusively superior to other methods. We recommend applying several different methods to the same data and comparing the resulting velocity models. The inversion of dispersion curves is essentially non-unique. It is important to understand that the results include uncertainty, and practitioners should try to incorporate other information in processing and interpretation to mitigate uncertainty to the extent possible. Applicability of the method and quality of data are very dependent on the site conditions, and we caution that the MAM is not applicable everywhere. We recommend initial processing of data immediately after data collection to verify that the raw data include sufficient coherent information. Unlike active seismic methods, the raw passive data of MAM look like noise, and it is difficult to evaluate data quality without spectral and/or coherency analysis. The calculation of phase velocity from ambient noise is not intuitive compared with other seismic methods. An understanding of fundamental theory and a sharing of experimental knowledge are important to ensure accuracy and reliability of the results of such investigations.

Appendices

Appendix 1

This appendix summarizes the essential features of the SPAC method. Let us start with two traces f(t) and g(t) obtained at two receives with separation Δx. For the sake of simplicity, waves are propagating parallel to a receiver array (Fig. 

Fig. 13

Spatial auto-correlation in 1D (a) and 2D (b). In the 1D wave propagation (a), waves propagate parallel to the receiver array, and cross-correlation between two receivers goes to cosine function. In the 2D array (b), cross-correlations are calculated between a receiver at the middle of the circle and receivers on the circle. The directional average of cosine function goes to Bessel function as shown in Fig. 14

13a). Two traces are transformed into frequency domain and written as F(ω) and G(ω). Cross-correlation of the two traces (CC_fg(ω)) can be defined as:

$${CC}_{fg}\left(\omega \right)=F\left(\omega \right)\cdot \overline{G\left(\omega \right)}={A}_{f}\left(\omega \right){A}_{g}\left(\omega \right)\cdot {\mathrm{exp}}^{i\Delta \varnothing \left(\omega \right)}$$

(6)

where A_f(ω) and A_g(ω) are amplitude of F(ω) and G(ω), respectively. Phase velocity c(ω) is directly related to the phase difference Δϕ(ω).

$$\Delta \varnothing \left(\omega \right)=\frac{\omega \cdot \Delta x}{c\left(\omega \right)}$$

(7)

Using Eq. (7), Eq. (6) can be rewritten as,

$${CC}_{fg}\left(\omega \right)=F\left(\omega \right)\cdot \overline{G\left(\omega \right)}={A}_{f}\left(\omega \right){A}_{g}\left(\omega \right)\cdot {\mathrm{exp}}^{i\left(\frac{\omega \cdot \Delta x}{c\left(\omega \right)}\right)}$$

(8)

Dividing the cross-correlation of the two traces (CC_fg(ω)) by autocorrelations (A_f(ω) and A_g(ω)) yields complex coherency of two traces COH_fg(ω).

$$CO{H}_{fg}\left(\omega \right)=\frac{C{C}_{fg}\left(\omega \right)}{\sqrt{{A}_{f}\left(\omega \right){A}_{g}\left(\omega \right)}}.$$

(9)

Using Eq. (8), Eq. (9) can be rewritten as,

$${\mathrm{COH}}_{fg}\left(\omega \right)=\frac{C{C}_{fg}\left(\omega \right)}{\sqrt{{A}_{f}\left(\omega \right){A}_{g}\left(\omega \right)}}={exp}^{i\left(\frac{\omega \bullet \Delta x}{c\left(\omega \right)}\right)},$$

(10)

and taking the real part of this gives

$$\mathrm{Re}\left({COH}_{fg}\left(\omega \right)\right)=\mathrm{cos}\left(\frac{\omega \cdot \Delta x}{c\left(\omega \right)}\right).$$

(11)

Equation (11) expresses the fact that the real part of the complex coherency for two traces takes the form of a cosine function. The SPAC coefficient for two-dimensional array data is defined as the directional average of complex coherencies Eq. (10), giving

$$SPAC\left(r,\omega \right)=\frac{1}{2\pi }{\int }_{\varphi =0}^{\varphi =2\pi }COH\left(r,\varphi ,\omega \right)d\varphi ,$$

(12)

where r is the separation of two receivers (or radius of a circle), and φ is the angular separation of the two sensors.

In essence, Eq. (12) calculates the complex coherencies for two sensors with separation r and angular separation φ and it averages these coherencies around a circle (Fig. 13b). The directional average of trigonometric functions reduces to a Bessel function as shown by Okada (2003):

$${J}_{0}\left(kr\right)=\frac{1}{2\pi }{\int }_{\varphi =0}^{\varphi =2\pi }\mathrm{cos}\left(kr,\varphi \right)d\varphi ,$$

(13)

where k is the wave number.

Figure 

Fig. 14

Comparison of cosine (cos(ω, φ)) functions corresponding to sensor pairs between center and A to G with orientation φ as shown in in Fig. 13 and the Bessel (J₀(ω)) function

14 compares cosine functions corresponding to sensor pairs between the center and receivers A to G with orientation φ as shown in Fig. 13, and Bessel functions. Using Eq. (13), the directional average of the right side of Eq. (12) can be written as.

$${J}_{0}\left(\frac{\omega \cdot r}{c\left(\omega \right)}\right)=\frac{1}{2\pi }{\int }_{\varphi =0}^{\varphi =2\pi }\mathrm{cos}\left(\frac{\omega \cdot r}{c\left(\omega \right)},\varphi \right)d\varphi ,$$

(14)

Finally, Eqs. (12) and (14) can be combined to get

$$\mathrm{Re}\left(SPAC\left(R,\omega \right)\right)={J}_{0}\left(\frac{\omega }{c\left(\omega \right)}r\right),$$

(15)

where c(ω) is the fundamental mode of phase velocity at the angular frequency ω.

When we process surface waves, seismic interferometry (SI) is essentially the same as spatial autocorrelation (SPAC). SPAC or coherency can be considered the real part of cross-correlation divided by autocorrelation in frequency domain. In other words, the inverse Fourier transform of SPAC corresponds to the cross-correlation in time domain. Figure 

Fig. 15

a Example of SPAC in the frequency domain calculated from microtremor measurement. b Phase velocity image in the frequency domain calculated from a. c SPAC in the time domain calculated from a by inverse Fourier transform. d Phase velocity image in the frequency domain calculated from c. (Modified from AH18)

15 shows an example of SPAC in the frequency domain (Fig. 15a) and time domain (Fig. 15c) calculated from microtremors (AH18 Section 10). Figure 15b shows the error between the coherencies and Bessel functions, and phase velocities that provide minimum error were defined as a dispersion curve. Frequency-domain coherencies (Fig. 15a) were transformed to time domain by the inverse Fourier transform (Fig. 15c). The time domain data are termed the SI or cross-correlation and resemble a shot gather. Further processing using phase shift and stack (Park et al. 1999) transforms it to a phase velocity versus frequency domain plot (Fig. 15d), and phase velocities that provide maximum stacked amplitude are defined as a dispersion curve. Figure 15a–d correspond to SPAC and SI, respectively, and it is evident that the dispersion curves obtained by two methods are almost identical.

Appendix 2

This appendix describes a traditional iterative non-linear least squares method in which the number and thickness of subsurface layers are fixed through the inversion process, and the layer specific V_S is the only unknown. As an alternative, both V_S and thickness of each layer are often estimated in the inversion. V_P and density are usually related to V_S linearly with empirical equations (e.g., Kitsunezaki et al. 1990; Ludwig et al. 1970) in each step of the iteration. The computation procedure is summarized as follows.

First, the set of Vs values in each layer can be written as a one-dimensional vector x, as.

$${x}^{T}=\left(V{s}_{1},V{s}_{2},\cdot \cdot \cdot \cdot ,V{s}_{M}\right),$$

(16)

where Vs₁,Vs₂, ・・・Vs_M are the V_S values for the 1st, 2nd, and Mth layers, respectively.

As an objective function, we choose the L2 norm between the observed and calculated phase velocity values from x, that is,

$${\sum }_{i}^{N}{\left({f}_{i}^{\mathrm{obs}}-{f}_{i}^{\mathrm{cal}}\left(V{s}_{1},V{s}_{2},\cdot \cdot \cdot \cdot V{s}_{M}\right)\right)}^{2}={\sum }_{i}^{N}{\left({f}_{i}^{\mathrm{obs}}-{f}_{i}^{\mathrm{cal}}\left(x\right)\right)}^{2}\to \mathrm{Minimum}$$

(17)

where N is the number of observed phase velocity data, f^obs are phase velocities obtained from observed waveform data and f^cal are theoretical phase velocities for the V_S model.

A theoretical phase velocity curve can be calculated by the matrix method (e.g., Saito and Kabasawa (1993). Now, a Jacobian matrix can be written as

$$a=\left(\begin{array}{c}\begin{array}{c}\begin{array}{ccc}\frac{\partial {f}_{1}^{cal}\left(x\right)}{{\partial Vs}_{1}}& \frac{\partial {f}_{1}^{cal}\left(x\right)}{{\partial Vs}_{2}}& \begin{array}{cc}\cdot & \frac{\partial {f}_{1}^{cal}\left(x\right)}{{\partial Vs}_{M}}\end{array}\end{array}\\ \begin{array}{ccc}\frac{\partial {f}_{2}^{cal}\left(x\right)}{{\partial Vs}_{1}}& \frac{\partial {f}_{2}^{cal}\left(x\right)}{{\partial Vs}_{1}}& \begin{array}{cc}\bullet & \frac{\partial {f}_{2}^{cal}\left(x\right)}{{\partial Vs}_{M}}\end{array}\end{array}\\ \begin{array}{ccc}\frac{\partial {f}_{3}^{cal}\left(x\right)}{{\partial Vs}_{1}}& \begin{array}{cc}\frac{\partial {f}_{3}^{cal}\left(x\right)}{{\partial Vs}_{2}}& \cdot \end{array}& \frac{\partial {f}_{3}^{cal}\left(x\right)}{{\partial Vs}_{M}}\end{array}\end{array}\\ \begin{array}{ccc}\cdot & \cdot & \begin{array}{cc}\cdot & \cdot \end{array}\end{array}\\ \begin{array}{ccc}\frac{\partial {f}_{N}^{cal}\left(x\right)}{{\partial Vs}_{1}}& \begin{array}{cc}\frac{\partial {f}_{N}^{cal}\left(x\right)}{{\partial Vs}_{2}}& \cdot \end{array}& \frac{\partial {f}_{N}^{cal}\left(x\right)}{{\partial Vs}_{M}}\end{array}\end{array}\right).$$

(18)

We can see that the unknown vector x is in the derivatives, and it makes the inversion non-linear. One of the consequences is that an initial guess model must be constructed, which must predict the data sufficiently well to obtain convergence. In actual calculation, elements of the Jacobian matrix a are usually calculated numerically using a finite-difference method (Xia et al. 1999). A residual between observed and theoretical phase velocities can be expressed as vector y given by:

$${y}^{T}=\left(\begin{array}{c}\begin{array}{c}{f}_{1}^{obs}-{f}_{i}^{cal}\left(x\right)\\ {f}_{2}^{obs}-{f}_{2}^{cal}\left(x\right)\\ {f}_{3}^{obs}-{f}_{3}^{cal}\left(x\right)\end{array}\\ \cdot \\ {f}_{N}^{obs}-{f}_{N}^{cal}\left(x\right)\end{array}\right)$$

(19)

The correction vector Δx to a given guess for f(x) can be calculated by the least squares method as follows:

$$\left({a}^{T}a+\varepsilon I\right)\Delta x={a}^{T}y,$$

(20)

where ε is a damping parameter stabilizing the inversion. In the lth iteration, a new estimated model x^l+1 is calculated as

$${x}^{l+1}={x}^{l}+\gamma \Delta x,$$

(21)

where γ is the constant specifying the step length. The method used to construct the initial model (x⁰) is described in Section 7.1.

Generally, the non-linear and non-unique nature of geophysical inverse problems requires the introduction of spatial regularization to converge to a physically meaningful result. In this case, it is common to regularize the variance in V_S between each layer as follows:

$$\left({a}^{T}a+\alpha {r}_{v}^{T}{r}_{v}+\varepsilon I\right)\Delta x={a}^{T}y,$$

(22)

where r_v is the difference between V_S of two successive layers. In matrix form, this equation can be expressed as

$${r}_{v}x=\left(\begin{array}{ccc}\begin{array}{c}1\\ 0\\ \begin{array}{c}0\\ \cdot \\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}-1\\ 1\end{array}\\ \begin{array}{c}0\\ \cdot \end{array}\\ 0\end{array}& \begin{array}{ccc}\begin{array}{c}0\\ \begin{array}{c}-1\\ 1\\ \cdot \end{array}\\ 0\end{array}& \begin{array}{c}\cdot \\ \begin{array}{c}\cdot \\ \cdot \\ \cdot \end{array}\\ \cdot \end{array}& \begin{array}{cc}\begin{array}{c}0\\ \begin{array}{c}0\\ 0\\ \cdot \end{array}\\ 1\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ \cdot \\ -1\end{array}\end{array}\end{array}\end{array}\right)\left(\begin{array}{c}\begin{array}{c}\begin{array}{c}{Vs}_{1}\\ {Vs}_{2}\end{array}\\ {Vs}_{3}\\ \cdot \end{array}\\ {Vs}_{M-1}\\ {Vs}_{M}\end{array}\right)=\left(\begin{array}{c}\begin{array}{c}0\\ 0\\ 0\end{array}\\ 0\\ 0\end{array}\right),$$

(23)

where α is the weight of regularization, and large α makes an inverted model smoother.