# Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays

Marc Wathelet, B Guillier, P Roux, C Cornou, M Ohrnberger

► To cite this version:

Marc Wathelet, B Guillier, P Roux, C Cornou, M Ohrnberger. Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays. Geophysical Journal International, 2018, 215 (1), pp.507-523. 10.1093/gji/ggy286. hal-04908196

HAL Id: hal-04908196

https://hal.science/hal-04908196v1

Submitted on 4 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Rayleigh wave three-component beamforming: signed ellipticity assessment from high-resolution frequency-wavenumber processing of ambient vibration arrays

M. Wathelet,<fnref n="1" /> B. Guillier,<fnref n="1" /> P. Roux,<fnref n="1" /> C. Cornou<fnref n="1" /> and M. Ohrnberger<fnref n="2" />

Accepted 2018 July 14. Received 2018 June 18; in original form 2018 April 05

## SUMMARY

The variation of Rayleigh ellipticity versus frequency is gaining popularity in site characterization. It becomes a necessary observable to complement dispersion curves when inverting shear wave velocity profiles. Various methods have been proposed so far to extract polarization from ambient vibrations recorded on a single three-component station or with an array of three-component sensors. If only absolute values were recovered 10 yr ago, new array-based techniques were recently proposed with enhanced efficiencies providing also the ellipticity sign. With array processing, higher-order modes are often detected even in the ellipticity domain. We suggest to explore the properties of a high-resolution beamforming where radial and vertical components are explicitly included. If *N* is the number of three-component sensors, $2N \times 2N$ cross-spectral density matrices are calculated for all presumed directions of propagation. They are built with *N* radial and *N* vertical channels. As a first approach, steering vectors are designed to fit with Rayleigh wave properties: the phase shift between radial and vertical components is either $-\Pi/2$ or $\Pi/2$. We show that neglecting the ellipticity tilt due to attenuation has only minor effects on the results. Additionally, we prove analytically that it is possible to retrieve the ellipticity value from the usual maximization of the high-resolution beam power. The method is tested on synthetic data sets and on experimental data. Both are reference sites already analysed by several authors. A detailed comparison with previous results on these cases is provided.

**Key words:** Fourier analysis; Time-series analysis; Site effects; Surface waves and free oscillations; Wave propagation.

## 1 INTRODUCTION

The dispersive property of surface waves has been used for decades to invert shear wave ($V_s$) profiles (Tokimatsu 1997). Over a wide range of frequency, ambient vibration arrays have proved their reliability for measuring dispersion curves either for Love and Rayleigh waves (Garofalo et al. 2016a). More recently, the ellipticity of Rayleigh waves is also gaining success with substantial improvements of the signal processing techniques. Hobiger et al. (2012) provide a complete and detailed overview of the methods developed so far to retrieve the experimental ellipticities, based on single three-component stations and arrays of three-component sensors.

For ambient vibrations, without any prior knowledge about the direction of propagation, a single station processing measures at best the absolute value of the ellipticity in a restricted frequency range, generally when the fundamental mode dominates the wavefield (Fäh et al. 2008; Hobiger et al. 2009). With a pair of three-component

sensors, Roueff et al. (2009) proposed a surface wave separation method when several polarized waves are recorded. With earthquake records on an array of three-component sensors, Boore & Toksöz (1969) already raised the interest of combining dispersion and ellipticity curve to resolve model ambiguities, even without considering the ellipticity sign. In a more systematic way, Hobiger et al. (2013) investigated the advantages and the limits of the inversion of the absolute value for a number of synthetic and experimental cases. Boaga et al. (2013) also highlighted the benefits of combining dispersion and ellipticity curves to correctly identify modes around osculation points.

Three-component array processing may improve the ellipticity resolution and bring new insights over the sign of the ellipticity and its variation versus frequency: prograde and retrograde transitions are critical observables that help to resolve model uncertainties (e.g. Scherbaum et al. 2003). Poggi & Fäh (2010) showed that the absolute value can be obtained from the ratio of the high-resolution

<footnote n="1">¹Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTERre, 38000 Grenoble, France. E-mail: Marc.Wathelet@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²Institute of Earth and Environmental Science, University of Potsdam, 14476 Potsdam, Germany</footnote>
beam power calculated on the radial and on the vertical compo-
nents. In favourable conditions, higher modes and their ellipticity
are also measured. MUSIQUE (Hobiger et al. 2012) and WaveDec
(Maranò et al. 2012, 2017) are currently the two competing methods
to calculate experimental ellipticity signs from ambient vibrations.
For MUSIQUE, slowness and direction of propagation are obtained
from a MUSIC (Schmidt 1986) analysis of a cross-spectral ma-
trix made by the sum of the one-component covariance matrices
 vertic, north, and east components). Horizontal signals are then
 projected on the transverse and the radial directions. If the radial en-
 ergy is higher than the transverse energy, a covariance matrix in the
 quaternion space is built from the vertical and the radial components
 transformed into a single quaternion signal. MUSIQUE delivers the
 phase difference between Rayleigh components and their ampli-
tude ratio. Maranò et al. (2012) suggested an alternative method
 not based on the cross-spectral density matrix. For WaveDec, the
 wavefield is modelled by a sum of waves each characterized by
 an amplitude, a phase and a wavenumber vector. Rayleigh waves
 have an additional parameter for their ellipticity. Parameters are
 estimated using an iterative process where the likelihood of the
 statistical model versus the three-component measurements is max-
 imized.

In the aforementioned references based on high-resolution beam-
forming (Capon 1969), cross-spectral matrices have always N × N
elements even if three-component wavefields are considered. We
propose to study the improvements that may be brought by 2N
× 2N cross-spectral matrices, made of radial and vertical compo-
nents under the assumption of Rayleigh plane waves. The idea of
combining all components into a single cross-spectral matrix was
first suggested by Wagner (1996). The steering matrices were con-
structed to extract polarization properties in the three directions.
Only phase shifts were considered. In addition, we are also check-
ing the feasibility of extracting the amplitude ratio between radial
and vertical components.

We first review the properties of conventional beamforming for
one to three components, before providing a high-resolution im-
plementation for three-component surface waves. For all of these
beamformers, analytical formulation is developed for a single plane
wave. A special attention is paid to the incoherent noise level that
may have a direct influence over quality of the ellipticity determi-
nation. The methodology is applied to synthetic wavefields and to
one experimental case.

2 METHODS

2.1 Conventional vertical beamforming (CVBF)

Assuming an array of *N* vertical sensors located on a horizontal
plane at $\vec{r}_i$ where *i* varies from 1 to *N* and a single harmonic plane
wave $S_{R_z}(\omega)$ ($\omega$ is the angular frequency, subscript $R_z$ stands for
the vertical component of a Rayleigh wave) which propagates with
wavenumber $\vec{k}_0$, the Fourier transforms of the signals recorded at
each sensor are defined in eq. (1). $X_Z$ is a $N \times 1$ column matrix.
Bold face characters are matrices.

$$
X_Z(\omega, \vec{k}_0) = S_{R_z}(\omega)q(\vec{k}_0) + \eta_Z(\omega), \quad (1)
$$

where $\eta_Z$ is the noise observed at each sensor which is incoherent
between any pair of sensors; $q(\vec{k})$ is a $N \times 1$ column matrix con-
taining the phase shifts for each sensor for a wave propagation at
wavenumber $\vec{k}$. It is defined by

$$
q(\vec{k}) = \exp(-j\vec{k} \cdot \vec{r}_i). \tag{2}
$$

Note that $q^*q = N$. Eq. (1) simply means that the same plane wave
is recorded at all sensors but with specific delays linked to the
magnitude and the sign of the projection of the propagation vector
on the position vector. The frequency wavenumber (f-k) power
spectral density $P(\omega, \vec{k})$ can be estimated with $\hat{P}(\omega, \vec{k})$ (Capon
1969). For the sake of simplicity, we drop the hat in this paper.

$$
P(\omega, \vec{k}) = \frac{1}{N^2} E^*(\vec{k}) F(\omega) E(\vec{k}), \quad (3)
$$

where

$$
F(\omega) = \frac{1}{M} \sum_{m=1}^{M} [X_{Z,m} X_{Z,m}^*]. \quad (4)
$$

$F$ is the cross-spectral density matrix estimated with the block
averaging method. $M$ is the number of non-overlapping blocks; $X_{Z,m}$
is the realization of $X_{R_z}$ for block $m$; $E(\vec{k}) = q(\vec{k})$ is the steering
matrix ($N \times 1$). If the noise $\eta_{R_z}$ is incoherent between any pair of
sensors and if $M$ is large enough, we have according to Capon (1969,
eq. 30),

$$
F(\omega) = |S_{R_z}|^2 (q_0 q_0^* + R I) = |S_{R_z}|^2 F_0, \quad (5)
$$

where $q_0 = q(\vec{k}_0)$; $R$ is the average ratio between the incoherent and coherent signal assumed to be the same for all sensors. Compared to Capon’s formulation, $F$ is not normalized to keep track of the original signal amplitudes along the whole process. With undistorted amplitudes, we show in the next sections that it is possible to extract the ellipticity information for three-component sensors directly from the maximization of an f-k power spectral density. Consequently, we assume that all sensors are perfectly calibrated. Combining eqs (3) and (5), $P$ reduces to

$$
\begin{align}
P(\omega, \vec{k}) &= \frac{|S_{R_z}|^2}{N^2} \left[ \mathbf{q}^* \mathbf{q}_0 \mathbf{q}_0^* \mathbf{q} + N R \right] \tag{6} \\
&= |S_{R_z}|^2 \left[ |B(\vec{k} - \vec{k}_0)|^2 + \frac{R}{N} \right], \nonumber
\end{align}
$$

where $|B(\vec{k})|^2$ is the beamforming array response pattern normalized
by the number of sensors $N$ and computed only from the array
geometry $\vec{r}_i$ (Capon 1969; Woods & Lintz 1973; Wathelet et al.
2008). It has a maximum value equal to unity at the origin, that
is, when $\vec{k}$ is equal to $\vec{k}_0$. However for a finite number of sensors,
$|B(\vec{k})|^2$ has also many side lobes whose amplitude may be close to
the central peak. If $\vec{k} = \vec{k}_0$, we have

$$
P(\omega, \vec{k}_0) = |S_{R_z}|^2 \left( 1 + \frac{R}{N} \right). \qquad (7)
$$

As a reference for the following discussion, Fig. 1 shows an example
of the response for an array made of 11 vertical sensors distributed
on a circle (100 m radius) and one at the centre of the circle [$k_{\min}$ =
0.024 rad m⁻¹ and $k_{\max}$ = 0.45 rad m⁻¹, see Wathelet et al. (2008)
for definitions]. The wavefield is composed of a single 10 Hz wave
travelling at 300 m s⁻¹ in the east direction. The beam pattern is
centred around $\vec{k}$ = [0.21, 0] with a radial symmetry in agreement
with the array layout.

2.2 Conventional radial beamforming (CRBF)

In this section we develop the radial projection for a conventional
beamforming. It follows the original ideas proposed by Fäh et al.
(2008) and Poggi & Fäh (2010) for a high-resolution beamforming.
The f-k power spectrum is also computed with N sensors but the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Vertical f-k power P obtained for a circular array made of 11 sensors on a circle with a radius of 100 m and one at the centre. A 10 Hz plane wave is travelling at 300 m s⁻¹ (grey circle) in the east direction (grey line). The small black circle indicates the wavenumber with the highest power. All points outside of a radius of 0.45 rad m⁻¹ are partially masked. It corresponds to the aliasing limit of the array (kmax).</figcaption></figure>

vertical component is replaced by a projection of the horizontal east and north components (usually aligned to the local magnetic north) on the direction of propagation (k̄).

We consider a wavefield made of a Rayleigh wave and a Love wave propagating in a 1-D layered medium. Rayleigh wave is distributed on the vertical component (Rz) and on the radial component (Rh) while the Love wave is only on the horizontal transverse component (L). The Rayleigh and Love waves travel at wavenumbers k̄₀ and k̄₁, respectively. Hence, we can define X<sub>Rh</sub> and X<sub>L</sub> as the Fourier transforms of the horizontal Rayleigh wave and the Love wave, respectively.

$$
\begin{align}
\mathbf{X}_{\mathrm{Rh}}(\omega, \vec{k}_0) &= S_{\mathrm{Rh}}(\omega)\mathbf{q}(\vec{k}_0) = S_{\mathrm{Rh}}(\omega)\mathbf{q}_0, \tag{8} \\
\mathbf{X}_{\mathrm{L}}(\omega, \vec{k}_1) &= S_{\mathrm{L}}(\omega)\mathbf{q}(\vec{k}_1) = S_{\mathrm{L}}(\omega)\mathbf{q}_1. \nonumber
\end{align}
$$

$S_{Rh}$ and $S_{Rz}$ are linked by the Rayleigh ellipticity $e_0$. According to Boore & Toksöz (1969), for perfectly elastic materials, Rayleigh ellipticity is an imaginary number either positive (prograde) or negative (retrograde). Inelasticity induces a tilt of the ellipticity with a small real part. For an extreme case ($Q_p = 10, Q_s = 4$ and a Poisson's ratio equal to 0.28), Boore & Toksöz (1969) found that the tilt is only 2 deg. From the equations provided by these authors, when Poisson's ratio tends to 0.5, the tilt vanishes. If we consider the general case where the ellipticity of the incoming wave is a complex number with a non-null real part, we have

$$
\frac{S_{\mathrm{Rh}}(\omega)}{S_{\mathrm{Rz}}(\omega)} = -e_0(\omega)(\tau_0(\omega) + j), \quad \tau_0(\omega) \in \mathbb{R}, e_0(\omega) \in \mathbb{R}, \quad (9)
$$

where $e_0(\omega)$ is the imaginary part of the ellipticity and $\tau_0(\omega)$ is the
ratio between the real and the imaginary parts. The dependency to $\omega$
is dropped further in the text for conciseness. For a tilt of 2 deg, $\tau_0$
$= 0.035$. Contrary to the usual definition found in Boore & Toksöz (1969), a minus sign is present in front of the ellipticity $e_0$ in our
formulation. $S_{Rh}$ is oriented in the direction of propagation and $S_{Rz}$
is pointing upwards to match with the usual convention followed by
the sensor manufacturers (e.g. Lennartz electronic GmbH), while a
vertical axis pointing downwards is generally adopted in the liter-
ature about theoretical dispersion curve computation (e.g. Dunkin
1965). It can be easily checked that a negative $e_0$ in eq. (9) produces
a retrograde elliptical particle motion.

For both experimental and numerically simulated wavefields, the
radial component is not obtained through a direct observation. In-
stead, it is mixed with transverse motion in horizontal components
that may have an arbitrary orientation compared to the Rayleigh
wave direction of propagation. We define the projections of the
Rayleigh and the Love waves on the east and the north components:

$$
\begin{align}
\mathbf{X}_E &= \cos(\theta_0)\mathbf{X}_{Rh} - \sin(\theta_1)\mathbf{X}_L + \eta_E \tag{10} \\
\mathbf{X}_N &= \sin(\theta_0)\mathbf{X}_{Rh} + \cos(\theta_1)\mathbf{X}_L + \eta_N,
\end{align}
$$

where η<sub>E</sub> and η<sub>N</sub> are some random noise observed at each sensor
on the horizontal components; θ<sub>0</sub> and θ<sub>1</sub> are the propagation angles
for the Rayleigh and the Love waves, respectively. The radial com-
ponents are estimated by projecting X<sub>E</sub> and X<sub>N</sub> on any arbitrary
direction θ = atan2(k<sub>y</sub>, k<sub>x</sub>).

$$
\begin{align}
\hat{\mathbf{X}}_{\mathrm{Rh}} &= \cos(\theta)\mathbf{X}_E + \sin(\theta)\mathbf{X}_N \nonumber \\
&= \mathbf{X}_{\mathrm{Rh}} \frac{\vec{k}_0 \cdot \vec{k}}{|\vec{k}_0||\vec{k}|} + \mathbf{X}_L \frac{(\vec{k}_1 \times \vec{k})_z}{|\vec{k}_1||\vec{k}|} + \frac{k_x}{|\vec{k}|} \eta_E + \frac{k_y}{|\vec{k}|} \eta_N \tag{11} \\
&= \alpha \mathbf{X}_{\mathrm{Rh}} + \alpha_L \mathbf{X}_L + \alpha_E \eta_E + \alpha_N \eta_N. \nonumber
\end{align}
$$

A hat is added to $\hat{\mathbf{X}}_{\mathrm{Rh}}$ because the quantity defined in eq. (11) is not the pure radial component of the Rayleigh wave ($\mathbf{X}_{\mathrm{Rh}}$) but only an estimation that contains also a Love term. For conciseness, we define $\alpha(\vec{k}_0, \vec{k})$ and $\alpha_L(\vec{k}_1, \vec{k})$ the projection functions for the Rayleigh and Love waves, respectively.

The cross-spectral density matrix is constructed in the same way as for a one-component beamforming, replacing $\mathbf{X}_{Rz}$ by $\hat{\mathbf{X}}_{\mathrm{Rh}}$. From the results developed in Appendix A, we have

$$
\begin{align}
F_{\text{Radial}}(\omega, \vec{k}) &= \frac{1}{M} \sum_{m=1}^{M} \hat{X}_{\text{Rh},m} \hat{X}_{\text{Rh},m}^* \tag{12} \\
&= \alpha^2 e_0^2 (1 + \tau_0^2) |S_{\text{Rz}}|^2 (F_0 + \beta F_1), \nonumber
\end{align}
$$

where

$$
\begin{align}
\beta(\vec{k}) &= \frac{\alpha_L(\vec{k})^2 |S_L|^2}{\alpha(\vec{k})^2 e_0^2 (1 + \tau_0^2) |S_{Rz}|^2}, \tag{13} \\
F_1 &= q_1 q_1^* + R I. \nonumber
\end{align}
$$

The conditions specified for eq. (5) also apply here. Additionally, the noise η<sub>E</sub> and η<sub>N</sub> must be uncorrelated for any pair of sensors including two components at the same location. The average ratios of incoherent over coherent signal at each sensor are the same for all components. F<sub>Radial</sub> may be estimated for a discrete number of directions (k̄) to save computation time (e.g. 72 directions over 2π). This approximation introduces very limited errors given that cos(2π/72) ≈ 0.99619. The resulting uncertainty is of the same order as the error on sensor orientation for which the maximum acceptable level is around 10 deg (Fäh et al. 2008).

The radial f-k spectral density becomes:

$$
\begin{align}
P_{\text{Radial}}(\omega, \vec{k}) &= \frac{\alpha^2 e_0^2 (1 + \tau_0^2) |S_{\text{Rz}}|^2}{N^2} \left[ \mathbf{q}^* \mathbf{q}_0 \mathbf{q}_0^* \mathbf{q} + N R + \beta (\mathbf{q}^* \mathbf{q}_1 \mathbf{q}_1^* \mathbf{q} + N R) \right] \tag{14} \\
&= \alpha^2 e_0^2 (1 + \tau_0^2) |S_{\text{Rz}}|^2 \left[ |B(\vec{k} - \vec{k}_0)|^2 + \beta |B(\vec{k} - \vec{k}_1)|^2 + \frac{R}{N} (1 + \beta) \right]. \nonumber
\end{align}
$$

If $\vec{k}_1$ is aligned with $\vec{k}_0$ ($\alpha_L = 0$ and $\beta = 0$), we have

$$
P_{\text{Radial}}(\omega, \vec{k}_0) = e_0^2 (1 + \tau_0^2) |S_{\text{Rz}}|^2 \left( 1 + \frac{R}{N} \right). \quad (15)
$$
## 2.3 Conventional Rayleigh three-component beamforming (CRTBF)

In a similar way as in eq. (3) we define a f-k power spectral density $P_R$ specific to Rayleigh waves. For a Rayleigh wavefield, 2N components ($\hat{X}_{Rh}$ and $X_{Rz}$) are available instead of just N like the common beamforming technique presented in Section 2.1. The definition of $P_R$ is given in eq. (16).

$$P_R(\omega, \vec{k}, e) = \frac{1}{N^2} E_R^*(\vec{k}, e) F_R(\omega, \vec{k}) E_R(\vec{k}, e), \quad (16)$$

where $e \in \mathbb{R}$ is the presumed Rayleigh ellipticity. The steering matrices do not include any ellipticity tilt due to attenuation, though the incoming wavefield has one ($\tau_0$). We show in the following that $\tau_0$ has almost no influence over the final results. We can define the steering matrix $E_R$ including ellipticity in two ways: either shift the radial components or the vertical components. We call $E_{Rh}$ and $E_{Rz}$, the steering matrices attached to the shift of radial and vertical components, respectively.

$$E_{Rh} = \begin{pmatrix} \frac{-j}{e}q \\ q \end{pmatrix}, \quad E_{Rz} = \begin{pmatrix} q \\ jeq \end{pmatrix} \quad (17)$$

The cross-spectral density matrix is constructed block-wise with the same order as $E_R$, that is, first the radial components, followed by the vertical ones.

$$F_R(\omega, \vec{k}) = \frac{1}{M} \sum_{m=1}^{M} \left[ \begin{pmatrix} \hat{X}_{Rh,m}(\vec{k}) \\ X_{Z,m} \end{pmatrix} \begin{pmatrix} \hat{X}_{Rh,m}^*(\vec{k}) & X_{Z,m}^* \end{pmatrix} \right]. \quad (18)$$

$\hat{X}_{Rh,m}$ and $X_{Z,m}$ are column matrices built from the Fourier transform of block $m$ for the radial and the vertical components, respectively. Combining eqs (1), (8), (9), (11) and (18), we have

$$F_R(\omega) = |S_{Rz}|^2 \times \begin{pmatrix} \alpha^2 e_0^2(1 + \tau_0^2)(F_0 + \beta F_1) - \alpha e_0(\tau_0 + j)q_0 q_0^* \\ -\alpha e_0(\tau_0 - j)q_0 q_0^* & F_0 \end{pmatrix}. \quad (19)$$

The conditions specified for eq. (5) also apply here. Additionally, the noise $\eta_{Rz}$, $\eta_E$ and $\eta_N$ must be uncorrelated for any pair of sensors including two components at the same location. The average ratios of incoherent over coherent signal at each sensor are the same for all components. Substituting $F_R$ and $E_R$ in eq. (16) by their respective values defined in eqs (17) and (19), we have

$$P_{Rz}(\omega, \vec{k}, e) = |S_{Rz}|^2 \left[ \begin{aligned} & ((\alpha e_0 + e)^2 + \alpha^2 e_0^2 \tau_0^2) |B(\vec{k} - \vec{k}_0)|^2 \\ & + \beta \alpha^2 e_0^2 (1 + \tau_0^2) |B(\vec{k} - \vec{k}_1)|^2 \\ & + \frac{R}{N} (\alpha^2 e_0^2 (1 + \tau_0^2) (1 + \beta) + e^2) \end{aligned} \right] \quad (20) \\ P_{Rh}(\omega, \vec{k}, e) = \frac{P_{Rz}}{e^2},$$

using $E_{Rz}$ and $E_{Rh}$ as steering matrices, respectively. The relation between $P_{Rh}$ and $P_{Rz}$ is true for any wavefield and not only for the particular case considered here. If the Love wave ($\vec{k}_1$) travels in the same direction as the Rayleigh wave ($\vec{k}_0$), for example, produced by a common source in the far field, $\beta$ vanishes for any $\vec{k}$ aligned with $\vec{k}_0$ and $P_R$ can be approximated by the first term in the neighbourhood of the direction pointed by $\vec{k}_0$ and an incoherent noise term. The Rayleigh f-k power is thus the product of the beam pattern, maximum at $\vec{k} = \vec{k}_0$, and a factor controlled by $e$ and $\vec{k}$ (through $\alpha$). $\alpha$ is maximum and is equal to 1 when $\vec{k}$ is in the right direction.

<figure><img src="image_4.png" /><figcaption>Figure 2. Conventional Rayleigh three-component beamforming (CRTBF) for the same circular array as in Fig. 1. A 10 Hz Rayleigh wave travels at 300 m s⁻¹ in the east direction (grey line). The small black circle indicates the wavenumber with the highest power. All points outside of a radius of 0.45 rad m⁻¹ are partially masked as in Figs 1(a) and (b) Rayleigh f-k power $P_R$ for $e = -1$ and $e = 1$, respectively. (c) Rayleigh f-k power $P_{R\pm}$. (d) Sign of the ellipticity: grey is positive and white is negative.</figcaption></figure>

At $\vec{k} = \vec{k}_0$ and neglecting all Love wave contributions, we have

$$P_{Rz}(\omega, \vec{k}_0, e) = |S_{Rz}|^2 \times \left[ (e_0 + e)^2 + \tau_0^2 e_0^2 + \frac{R}{N} (e_0^2(1 + \tau_0^2) + e^2) \right]. \quad (21)$$

$P_R$ is maximum if the sign of $e$ is the same as the sign of the true ellipticity $e_0$, even if $R$ is large, if there are only a few sensors and/or if the real part of the true ellipticity ($\tau_0 e_0$) cannot be neglected. It is not possible to retrieve the ellipticity value when maximizing $P_R$ with any of the two proposed steering matrices $E_{Rh}$ and $E_{Rz}$. However, under the restriction that $e = \pm 1$, looking for the maximum of $P_R$ provides both $k_0$ and the sign of $e_0$. The two options for the steering matrix are equivalent if $e$ is forced to be $\pm 1$. To display the values of $P_R$, we may use $P_{R\pm}$ defined as:

$$P_{R\pm}(\omega, \vec{k}) = \max_{e=\pm 1} P_R(\omega, \vec{k}, e). \quad (22)$$

Fig. 2 shows $P_R$, $P_{R\pm}$, and the sign of the ellipticity $e$ for the same array as in Section 2.1 with a wavefield composed of a pure 10 Hz Rayleigh wave ($e_0 = 1$) travelling in the east direction at 300 m s⁻¹. Comparing Figs 1 and 2(c), the same pattern is observed except for side lobes that are not in the vicinity of the direction of propagation ($\vec{k}_0$). They are vanishing around the north axis as already observed from eq. (20). The main peak being in the middle of a large area with all positive ellipticities (Fig. 2d), the extraction of the ellipticity sign is robust against any small error of the velocity estimation.

Let us now consider a more realistic case including a Love wave travelling in the same direction as the Rayleigh wave but
at 320 m s⁻¹. CRBF and CRTBF are shown in Figs 3(a) and (b). The patterns look similar but the relative amplitudes of peaks are different. The location of the highest peak is close to the true wavenumber for CRTBF while it is erroneous for CRBF. The peak with the true wavenumber is much lower than the highest side lobes influenced by the dominant Love wave. The influence of the Love wave can be observed by comparing Figs 2(c) and 3(b) (peaks around north axis). The ellipticity sign pattern is sharper for the case without any Love wave (Figs 2d and 3c), but the sign remains correct in the vicinity of the beamforming peak when a Love wave is present. CRTBF is performing slightly better than CRBF. If we note that $P_{\text{radial}}$ is equal to $P_{Rz}$ at $e = 0$, the ellipticity $e$ adds a little more weight to the Rayleigh wave relative to the incoherent noise term ($R/N$) in eq. (20) compared to a radial beamforming (RBF). This advantage disappears if the amplitude of Love wave increases to 3 or more while Rayleigh amplitude remains at 1.

In Fig. 4, still with the same circular array, we check the capabilities of retrieving the correct properties of many waves travelling with similar velocities and whose azimuth are equally distributed from 90 to 140 deg. Love waves are travelling at 240 m s⁻¹, Rayleigh waves at 200 and 220 m s⁻¹ with ellipticities of -0.5 and 2, respectively. All sources have the same amplitude. The ratio of incoherent noise over coherent waves is fixed to 12 per cent. 48 blocks are averaged in all cases to compute the cross-spectral density matrix. For all cases the image appears relatively blurred. For CVBF (Fig. 4a), the two Rayleigh modes are visible but all incident plane waves are not fully recovered and aliasing side lobes have almost similar amplitudes even in the restricted area shown in Fig. 4. Many side lobes exist outside the visible area with the same intensity as the maxima in the visible area. CRBF (Fig. 4b) is much clearer with respect to side lobes even including those outside of the visible area, but the second Rayleigh mode (at 200 m s⁻¹) is not visible at all. Only the mode with the highest absolute value of ellipticity is evidenced. Positive CRTBF [$P_R$ with $e = 1$, Fig. 4(c)] looks very similar to CRBF, dominated by the Rayleigh wave with the strongest energy on the horizontal component (220 m s⁻¹ and $e = 2$). The side lobes outside of the visible area are smaller than for CRBF. If the ellipticity sign is forced to -1, Fig. 4(e) shows that the Rayleigh mode with the lowest absolute ellipticity (at 200 m s⁻¹) can still be extracted with a prior information about the direction. Strong side lobes are observed in the opposite direction, outside of the visible area. The ellipticity sign (Fig. 4e) confirms that the positive ellipticity dominates the f-k spectrum. Nevertheless negative patches are visible for all waves at 200 m s⁻¹. The same case will be processed with high-resolution techniques in the next sections.

In this section, we show that, even for a conventional beamforming, combining vertical and radial components inside the beamforming itself is a good way to retrieve the sign of the ellipticity. The assessment of phase slownesses is of the same quality as a usual one-component processing (vertical or radial). The CPU cost of such procedure remains reasonable. For instance, to process 1 hr of ambient vibrations with 14 sensors, it takes approximately 33 s for CVBF and about 100 s for CRTBF, keeping exactly the same parameters and on the same 16-core machine with similar CPU loads. The conventional beamforming used so far in this paper suffers strong limitations. Among others, we cannot determine the absolute value of ellipticity. Though it can be circumvented by combining our results with the absolute ellipticity values obtained with the method proposed by Poggi & Fäh (2010), we show in the next section that a high-resolution algorithm used with a combination

of vertical and radial components offers very interesting properties, solving both the sign and the absolute value issues for ellipticity.

## 2.4 High-resolution vertical beamforming (VBF)

According to Capon (1969) the high-resolution f-k power spectral density is

$$P'(\omega, \vec{k}) = \left[ E^*(\vec{k}) [F(\omega)]^{-1} E(\vec{k}) \right]^{-1}, \quad (23)$$

where $F$ is defined by eq. (4). The number $M$ of non-overlapping blocks used to compute $F$ must be larger than $N$, the number of sensors in the array.

For a wavefield composed of a plane wave propagating at $\vec{k}_0$, $F$ reduces to $F_0$ (eq. 5). The inverse of $F_0$ is based on Capon (1969, eq. 31).

$$F_0^{-1} = \frac{1}{R} \left[ I - \frac{q_0 q_0^*}{N + R} \right]. \quad (24)$$

Hence,

$$P'(\omega, \vec{k}) = |S_{Rz}|^2 \frac{R}{N} \left[ 1 - \frac{|B(\vec{k} - \vec{k}_0)|^2}{1 + \frac{R}{N}} \right]^{-1}, \quad (25)$$

which is equivalent to eq. (32) in Capon (1969). At $\vec{k} = \vec{k}_0$,

$$P'(\omega, \vec{k}_0) = |S_{Rz}|^2 \left(1 + \frac{R}{N}\right) = P(\omega, \vec{k}_0). \quad (26)$$

## 2.5 High-resolution radial beamforming (RBF)

As in Poggi & Fäh (2010), a high-resolution f-k power can be defined by injecting the definition of $F_{\text{Radial}}$ (eq. 12) in eq. (23). For the sake of simplicity, we do not consider any Love waves in the following analytical developments, but numerical tests will show that it has almost no influence over the determination of $\vec{k}$ in accordance with the conclusions drawn by Capon (1969) when two waves are present in the wavefield.

$$P'_{\text{Radial}}(\omega, \vec{k}) = \alpha^2 e_0^2 (1 + \tau_0^2) |S_{Rz}|^2 \frac{R}{N} \left[ 1 - \frac{|B(\vec{k} - \vec{k}_0)|^2}{1 + \frac{R}{N}} \right]^{-1}. \quad (27)$$

The absolute value of ellipticity can be obtained by dividing $P'_{\text{Radial}}$ by $P'$ (Poggi & Fäh 2010). From eq. (27), it comes that

$$|e| = \sqrt{\frac{P'_{\text{Radial}}}{P'}} = |\alpha(\vec{k})e_0| \sqrt{1 + \tau_0^2}. \quad (28)$$

For a 2 deg tilt, the square root equals to 1.0006, which may introduce a very small distortion. The obtained ellipticity is correct if $\vec{k}$ is in the same direction as $\vec{k}_0$ ($\alpha \approx 1$), even if the absolute value of the velocity estimation is false. The ellipticity estimator is relatively robust to small errors on the direction of propagation, for example, for an angle of 20 deg between $\vec{k}$ and $\vec{k}_0$, the error on the ellipticity is limited to about 5 per cent. In the method proposed by Poggi & Fäh (2010), $\vec{k}$ is retrieved independently from either the vertical or the radial component, and both results are generally plotted together. The quality of the ellipticity determination depends only upon the efficiency of these beamformers to properly decompose the wavefield and to extract the correct direction of propagation.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** CRBF and CRTBF for a wavefield composed of Love (320 m s⁻¹) and Rayleigh (300 m s⁻¹) waves with amplitudes 2 and 1, respectively. The array is the same circle as in Fig. 1. Phase differences between Love and radial Rayleigh wave are random across the 20 blocks used to build the cross-spectral density matrix. (a) CRBF. (b) CRTBF. (c) The ellipticity sign deduced from CRTBF (grey is positive, white is negative).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Comparison of several conventional beamforming techniques: (a) CVBF, (b) CRBF, (c) CRTBF with positive ellipticity ($P_R$ with $e = 1$), (d) CRTBF with negative ellipticity ($P_R$ with $e = -1$) and (e) ellipticity sign from $P_{R±}$ (grey is positive, white is negative). The level of incoherent noise is 12 per cent ($R/N = 0.01$). 48 blocks are used in all cases to evaluate the cross-correlation matrices. The wavefield is composed of six Rayleigh waves at 220 m s⁻¹ and $e = 2, 6$; Rayleigh waves at 200 m s⁻¹ and $e = -0.5$ and 6 Love waves at 240 m s⁻¹. All waves have a unitary intensity. The small subplots in panels (a) to (d) show the full wavenumber space between −0.33 and 0.33 rad m⁻¹. The colour scales are normalized in each plot to the maximum value of the f-k power spectrum in the visible area. The straight lines indicate the direction of propagation that vary from 90 to 140 deg. The circles are located at the theoretical velocities: 240, 220 and 200 m s⁻¹ (from left to right in the main plots).</figcaption></figure>
## 2.6 High-resolution Rayleigh three-component beamforming (RTBF)

As detailed in the previous section, the *N* radial and the *N* vertical components can be considered as part of a super array made of 2*N* single-component sensors. We propose to check the properties of a high-resolution Rayleigh beam power computed with 2*N* signals instead of just *N*. Following the high-resolution beam power defined by Capon (1969, eq. 18), we have

$$P'_{R}(\omega, \vec{k}, e) = \left[ \mathbf{E}'_{R}^{*}(\vec{k}, e) [\mathbf{F}_{R}(\omega)]^{-1} \mathbf{E}'_{R}(\vec{k}, e) \right]^{-1} \quad (29)$$

$\mathbf{F}_R$ is estimated by the means of eq. (19) with *M* non-overlapping blocks, *M* being larger than 2*N*. We can define the steering matrices $\mathbf{E}'_R$ including ellipticity in two ways like in the previous section. We call $\mathbf{E}'_{Rh}$ and $\mathbf{E}'_{Rz}$, the steering matrices attached to the shift of radial and vertical components, respectively.

$$\mathbf{E}'_{Rh} = \begin{pmatrix} -je\mathbf{q} \\ \mathbf{q} \end{pmatrix}, \quad \mathbf{E}'_{Rz} = \begin{pmatrix} \mathbf{q} \\ j/e \mathbf{q} \end{pmatrix}. \quad (30)$$

These matrices differ from those defined in eq. (17) because the high-resolution beam power is computed with the inverse of the cross-spectral density matrix. At $e = e_0$, the steering matrices in eq. (17) are equalizing the amplitudes of the vertical and the radial components of the cross-spectral density matrix. Intuitively, the radial components in a high-resolution beam power must be multiplied by the presumed ellipticity *e* instead of being divided to achieve a similar equalization. The proof is provided below.

For a wavefield composed of a Rayleigh plane wave propagating at $\vec{k}_0$ and with ellipticity $e_0$, we now show that $P'_R(\omega, \vec{k}, e)$ is maximum when $\vec{k} = \vec{k}_0$ and $e = e_0$. As in Section 2.5 no Love wave is considered. We develop a similar approach as Capon (1969, eq. 26) but generalized to three-component sensors and without $\mathbf{F}_R$ normalization. According to Henderson & Searle (1981, eq. 8), the block-wise inverse of $\mathbf{F}_R$ can be written as

$$\begin{aligned} F_R^{-1} &= \frac{1}{|S_{Rz}|^2} \begin{pmatrix} \alpha^2 e_0^2 (1 + \tau_0^2) F_0 & -\alpha e_0 (\tau_0 + j) q_0 q_0^* \\ -\alpha e_0 (\tau_0 - j) q_0 q_0^* & F_0 \end{pmatrix}^{-1} \\ &= \frac{1}{\alpha^2 e_0^2 (1 + \tau_0^2) |S_{Rz}|^2} \times \begin{pmatrix} F_0^{-1} + F_0^{-1} q_0 q_0^* A^{-1} q_0 q_0^* F_0^{-1} & \alpha e_0 (\tau_0 + j) F_0^{-1} q_0 q_0^* A^{-1} \\ \alpha e_0 (\tau_0 - j) A^{-1} q_0 q_0^* F_0^{-1} & \alpha^2 e_0^2 (1 + \tau_0^2) A^{-1} \end{pmatrix}, \end{aligned} \quad (31)$$

where $\mathbf{A} = \mathbf{F}_0 - q_0 q_0^* F_0^{-1} q_0 q_0^*$ and the inverse of $\mathbf{A}$ can be simplified to

$$A^{-1} = \frac{1}{R} \left[ \mathbf{I} - \frac{q_0 q_0^*}{2N + R} \right]. \quad (32)$$

Hence,

$$F_R^{-1} = \frac{1}{R \alpha^2 e_0^2 (1 + \tau_0^2) |S_{Rz}|^2} \times \begin{pmatrix} \left[ \mathbf{I} - \frac{\mathbf{q}_0 \mathbf{q}_0^*}{2N + R} \right] & \alpha e_0 (\tau_0 + j) \frac{\mathbf{q}_0 \mathbf{q}_0^*}{2N + R} \\ \alpha e_0 (\tau_0 - j) \frac{\mathbf{q}_0 \mathbf{q}_0^*}{2N + R} & \alpha^2 e_0^2 (1 + \tau_0^2) \left[ \mathbf{I} - \frac{\mathbf{q}_0 \mathbf{q}_0^*}{2N + R} \right] \end{pmatrix}. \quad (33)$$

Using $\mathbf{E}'_{Rh}$ in eq. (29), it comes that

$$P'_{Rh}(\omega, \vec{k}, e) = \frac{\alpha^2(1 + \tau_0^2) |S_{Rz}|^2}{N} \frac{\left( \frac{e_0^2}{e_0^2 + \alpha^2(1 + \tau_0^2)} - \left( \frac{e_0}{e_0^2 + \alpha^2(1 + \tau_0^2)} \right)^2 \frac{|B(\vec{k} - \vec{k}_0)|^2}{2 + \frac{R}{N}} \right)}{(\frac{e_0^2}{e_0^2} + \alpha^2(1 + \tau_0^2)) - \left( \frac{e_0}{e_0^2 + \alpha^2(1 + \tau_0^2)} \right)^2 \frac{|B(\vec{k} - \vec{k}_0)|^2}{2 + \frac{R}{N}}}. \quad (34)$$

In a similar way, using $E_{Rz}$ in eq. (29), we obtain

$$P'_{Rz}(\omega, \vec{k}, e) = e^2 P'_{Rh}(\omega, \vec{k}, e). \quad (35)$$

As in Section 2.3 the two steering matrices $E_{Rh}$ and $E_{Rz}$ lead to similar f-k power spectral densities for any wavefield. Hence, it is not necessary to compute the two steering matrices explicitly. $P'_{Rz}$ can be computed by multiplying $P'_{Rh}$ by the square of the presumed ellipticity.

The roots of the derivatives of $P'_{Rh}$ and $P'_{Rz}$ versus *e* correspond to maxima are found at

$$\begin{aligned} \hat{e}_h &= e_0 \alpha \frac{|B(\vec{k} - \vec{k}_0)|^2}{2 + \frac{R}{N} - |B(\vec{k} - \vec{k}_0)|^2} \\ \hat{e}_z &= e_0 \alpha \frac{(1 + \tau_0^2)[2 + \frac{R}{N} - |B(\vec{k} - \vec{k}_0)|^2]}{|B(\vec{k} - \vec{k}_0)|^2}, \end{aligned} \quad (36)$$

where $\hat{e}_h$ and $\hat{e}_z$ are the values of ellipticity inferred from the maximum of $P'_{Rh}$ and $P'_{Rz}$, respectively. When $\vec{k} = \vec{k}_0$ ($\alpha = 1$), $P'_{Rh}$ becomes

$$P'_{Rh}(\omega, \vec{k}_0, e) = |S_{Rz}|^2 \frac{R}{N} \left(2 + \frac{R}{N}\right) \times \frac{1 + \tau_0^2}{\left(\frac{e}{e_0} - 1\right)^2 + \tau_0^2 + \frac{R}{N}\left(\frac{e_0^2}{e_0^2} + 1 + \tau_0^2\right)}. \quad (37)$$

At $\vec{k} = \vec{k}_0$ and $e = e_0$, neglecting $\tau_0^2$, we have

$$P'_{Rh}(\omega, \vec{k}_0, e_0) \approx |S_{Rz}|^2 (1 + \frac{R}{2N}), \quad (38)$$

which is a clear reduction of the incoherent noise term by a factor 2 compared to $P'$ (eq. 26). The variation of $P'_R$ versus the presumed ellipticity *e* is shown in Fig. 5 for two true ellipticities: -0.5 (a) and 2 (b). All the continuous lines (low noise and/or high number of sensors) have a sharp and well-defined peak at the true ellipticity value. At $\vec{k} = \vec{k}_0$, the estimated ellipticities $\hat{e}_h$ and $\hat{e}_z$ are

$$\begin{aligned} \hat{e}_h &= \frac{e_0}{1 + \frac{R}{N}} \\ \hat{e}_z &= e_0 \left(1 + \frac{R}{N}\right) \left(1 + \tau_0^2\right). \end{aligned} \quad (39)$$

The position of the maximum of $P'_{Rh}$ is not influenced by the real part ($\tau_0$) of the true ellipticity. As shown in Fig. 5, $P'_{Rh}$ and $P'_{Rz}$ are almost not affected by an increase of the tilt from 0 to 2 deg (Boore & Toksöz 1969). In Fig. 5 and in eq. (39) the maxima of $P'_{Rh}$ (green) and $P'_{Rz}$ (orange) are moving in opposite directions when the noise increases, providing a biased estimation of $e_0$. The sharpness of the peak is also lost. Therefore the ratio between the results of two steering matrices is a marker of the noise level and we can define an estimated noise level by $\hat{R}$, neglecting any ellipticity tilt.

$$\hat{R} = N \left[ \sqrt{\frac{\hat{e}_z}{\hat{e}_h}} - 1 \right] \quad (40)$$

For experimental cases and noisy data, the true velocity may not be accurately resolved. In this case $\vec{k}$ is slightly different from $\vec{k}_0$ ($|B|^2 < 1$) in eq. (34). Let us consider the case $|B|^2 = 0.5$. For a direction of $\vec{k}$ in the neighbourhood of the direction of $\vec{k}_0$ so that $\alpha \approx 1$, the retrieved ellipticities becomes

$$\begin{aligned} \hat{e}_h &= e_0 \frac{1}{3 + 2 \frac{R}{N}} \\ \hat{e}_z &= e_0 (3 + 2 \frac{R}{N}) (1 + \tau_0^2). \end{aligned} \quad (41)$$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** $P'_{Rh}(\vec{k}_0, e)$ and $P'_{Rz}(\vec{k}_0, e)$ for (a) $e_0 = -0.5$ and (b) $e_0 = 2$. $P'_{Rh}$ is shown with blue ($\tau_0 = 0$) and green ($\tau_0 = 0.035$, a 2 deg tilt) lines. $P'_{Rz}$ is shown with red ($\tau_0 = 0$) and orange ($\tau_0 = 0.035$) lines. Continuous lines have a ratio $R/N = 0.01$ and dashed lines have a ratio $R/N = 0.1$. Blue and red curves are slightly visible only around the peaks, elsewhere they are almost the same as the green and the orange curves, respectively.</figcaption></figure>

$P'_{Rh}$ and $P'_{Rz}$ are no longer maximum for $e = e_0$ even for a very small amount of incoherent noise. As the distance between $\vec{k}$ and $\vec{k}_0$ increases, $|B|^2$ tends to 0 and $\hat{e}_h$ decreases to 0 while $\hat{e}_z$ tends to $\pm\infty$. Thus an accurate determination of $\vec{k}$ and a low level of incoherent noise are required to obtain an accurate ellipticity, which is likely to be impossible for real data sets. Alternatively, we propose to estimate the quantity $P'_{Rs}$ defined by the product of $P'_{Rh}$ and $P'_{Rz}$.

$$P'_{Rs}(\omega, \vec{k}, e) = P'_{Rh} P'_{Rz} = e^2 (P'_{Rh})^2, \qquad (42)$$

which has two maxima at the presumed ellipticities (see Appendix B)

$$\hat{e}_s = \pm|e_0||\alpha|\sqrt{1 + \tau_0^2}. \qquad (43)$$

This is the same formula as in eq. (28) for the ellipticity calculated from the ratio of the radial and vertical f-k power spectral densities. The ellipticity estimation is no longer influenced by the amount of incoherent noise or the number of sensors. As noted earlier when commenting about eq. (28), the alignment of $\vec{k}$ and $\vec{k}_0$ is the most important factor controlling the accuracy of the ellipticity measurement. Contrary to the method proposed by Poggi & Fäh (2010) where only the absolute value can be retrieved, we have here two maxima of the f-k power spectral density that may have distinct amplitudes. From eqs (37) and (42), we have at $\vec{k} = \vec{k}_0$

$$P'_{Rs}(\omega, \vec{k}_0, e) = |S_{Rz}|^4 \frac{R^2}{N^2} (2 + \frac{R}{N})^2 \times \frac{(1 + \tau_0^2)^2 e^2}{[(\frac{e}{e_0} - 1)^2 + \tau_0^2 + \frac{R}{N}(\frac{e_0^2}{e_0^2} + 1 + \tau_0^2)]^2}. \qquad (44)$$

At $|B|^2 = 0.5$ and assuming that $\vec{k}$ is still sufficiently aligned with $\vec{k}_0$ so that $\alpha \approx 1$,

$$P'_{Rs}(\omega, |B|^2 = 0.5, e) = |S_{Rz}|^4 \frac{R^2}{N^2} (2 + \frac{R}{N})^2 \times \frac{4(1 + \tau_0^2)^2 e^2}{[(\frac{e}{e_0} - 1)^2 + \tau_0^2 + (2 + 2\frac{R}{N})(\frac{e_0^2}{e_0^2} + 1 + \tau_0^2)]^2}. \qquad (45)$$

In both eqs (44) and (45), the sign of $e$ only influences the first term of the denominator, which is null for $e = e_0$ and positive when $e = -e_0$. Hence the maximum at $e = e_0$ is always higher than the

maximum at $e = -e_0$. The difference tends to zero when the noise level increases. For $|B|^2 = 0.5$, the noise term is larger than at $\vec{k} = \vec{k}_0$. An example of variations of $P'_{Rs}$ versus $e$ is given in Fig. 6 which confirms that the ellipticity value and its sign are robust against the amount of incoherent noise and the inaccuracy of $\vec{k}$. If equivalent maxima for the positive and the negative ellipticity are observed, it is a clear indication that the velocity measure is not accurate or that the noise level is too large.

In Fig. 7, with the same circular array and the same wavefield as in Fig. 4, we check the capabilities of retrieving the correct properties from high-resolution techniques with multiple sources. In Figs 7(a) and (b), the f-k power is a two-parameter function while a three-parameter function must be shown in Fig. 7(c). Hence, values of $P'_{Rxy}(\vec{k})$ defined by eq. (46) are displayed instead of a complex 3-D plot. In Fig. 6, $P'_{Rs}(\vec{k}_0)$ has always two maxima, one at the true ellipticity and the other at the opposite ellipticity. Keeping only the maximum value for a fixed $k_x$ and $k_y$ has the advantage of focusing only on the highest of the two peaks that is supposed to appear at the true ellipticity. These conclusions were drawn from the single source wavefield but we can extrapolate them to multiple sources in the neighbourhood of the true wavenumbers when the array properly separates the propagating waves.

$$P'_{Rxy}^{S}(\vec{k}) = \max_{e} P'_{Rs}(\vec{k}, e). \qquad (46)$$

Compared to conventional techniques, the peaks of the f-k power spectra are better defined and clearly higher than aliasing side lobes. For the vertical beamforming (VBF, Fig. 7a), the two Rayleigh modes are visible but many aliasing patterns are also present outside of the visible area. For the RBF (Fig. 7b) all peaks of the first Rayleigh mode (at 220 m s⁻¹) are clearly identified, but the second Rayleigh mode (at 200 m s⁻¹) is not visible at all. Only the mode with the highest ellipticity absolute value is evidenced as in Fig. 4(b). The Rayleigh beamforming (Fig. 7c) correctly identifies all waves for the two Rayleigh modes with a low level for all side lobes, in a better way than other techniques. The absolute values of the ellipticity are plotted in Figs 7(d) and (e), comparing results from Poggi's method and from RTBF, respectively.

Highest peaks are automatically extracted from the results in Figs 7(a)–(c) considering the complete wavenumber space. We start with a coarse grid search in $k_x$ and $k_y$ with the appropriate powers (for Figs 7a and b) or the function defined in eq. (46) (for Fig. 7c),
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Variations of $P'_s$ versus the presumed ellipticity $e$: (a) for $e_0 = -0.5$ and (b) for $e_0 = 2$. Plain and dashed lines have a ratio $R/N$ equal to 0.01 and 0.1, respectively. Black and grey lines are computed at $|B|^2$ equal to 1 ($\bar{k} = \bar{k}_0$) and 0.5, respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of several high-resolution beamforming techniques: (a) VBF, (b) RBF and (c) RTBF. (d) Ellipticities calculated with the method proposed by Poggi & Fäh (2010). (e) Absolute values of the ellipticities calculated with RTBF. The same wavefield as in Fig. 4 is used. The level of incoherent noise is 12 per cent ($R/N = 0.01$). 48 blocks are used in all cases to evaluate the cross-correlation matrices. The small subplots in panels (a)–(c) show the full wavenumber space between $-0.33$ and $0.33$ rad m$^{-1}$. The colour scales are normalized in each plot to the maximum value of the f–k power spectrum in the visible area. The straight lines indicate the direction of propagation that vary from 90 to 140 deg. The circles are located at the theoretical velocities: 240, 220 and 200 m s$^{-1}$ (from left to right in the main plots).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Histograms of velocity values picked in Fig. 7. (a) VBF, (b) RBF and (c) RTBF.</figcaption></figure>

<table><caption>**Table 1.** Properties of the power peaks automatically picked in Fig. 7(c) and part of the two main spikes in Fig. 8(c). Directions are counted clockwise from north (y axis).</caption><thead><tr><th>Power</th><th>Velocity (m s<sup>-1</sup>)</th><th>Direction (deg)</th><th>Ellipticity</th><th>Poggi Ell.</th></tr></thead><tbody><tr><td>3.44</td><td>220</td><td>90</td><td>1.88</td><td>1.87</td></tr><tr><td>2.73</td><td>220</td><td>100</td><td>1.90</td><td>1.99</td></tr><tr><td>3.86</td><td>220</td><td>110</td><td>1.92</td><td>2.08</td></tr><tr><td>2.75</td><td>220</td><td>120</td><td>1.94</td><td>1.98</td></tr><tr><td>3.19</td><td>220</td><td>130</td><td>1.92</td><td>2.01</td></tr><tr><td>3.37</td><td>220</td><td>140</td><td>1.87</td><td>1.68</td></tr><tr><td>0.26</td><td>200</td><td>90</td><td>-0.72</td><td>0.71</td></tr><tr><td>0.36</td><td>200</td><td>100</td><td>-0.73</td><td>0.76</td></tr><tr><td>0.32</td><td>200</td><td>110</td><td>-0.76</td><td>0.81</td></tr><tr><td>0.32</td><td>200</td><td>120</td><td>-0.72</td><td>0.77</td></tr><tr><td>0.31</td><td>200</td><td>130</td><td>-0.73</td><td>0.78</td></tr><tr><td>0.29</td><td>200</td><td>140</td><td>-0.69</td><td>0.75</td></tr><tr><td>0.37</td><td>201</td><td>387</td><td>-1.74</td><td></td></tr><tr><td>0.29</td><td>202</td><td>302</td><td>-1.74</td><td></td></tr></tbody></table>

itself calculated with a coarse sampling of arctan *e* from −Π/2 to Π/2 with a constant step of Π/20. Then, we refine the peak location in the 2-D ($k_x$ and $k_y$) or 3-D space ($k_x$, $k_y$ and arctan *e*) down to a relative precision of 10⁻⁵. Only the peaks with an amplitude above 5 per cent of the maximum amplitude are selected. Histograms of the corresponding velocity values are shown in Fig. 8. For RTBF two clear spikes are visible at 200 and 220 m s⁻¹ while no particular velocity can be identified from the vertical and radial histograms. All peaks which are part of the two spikes are selected and summarized in Table 1. All 12 velocities and their directions are retrieved with an accuracy of less than 1 m s⁻¹ and 1 deg, respectively. Two aliasing peaks are also extracted with matching velocities (last two lines in Table 1). RTBF correctly extracts the sign of the ellipticity for all 12 individual Rayleigh waves. The ellipticity of the waves at 200 m s⁻¹ (*e* = −0.5) and 220 m s⁻¹ (*e* = 2) are retrieved with a systematic relative error of about 50 and 5 per cent, respectively. The uncertainty is higher for the waves with the ellipticity having the lowest absolute value. Assuming that velocities and directions are perfectly determined, the ellipticities obtained by the method

proposed by Poggi & Fäh (2010) are compared with those obtained from RTBF in Table 1. In Fig. 7(e), the spots at 220 m s⁻¹ are larger than those found in Fig. 7(d). The ellipticities deduced from Fig. 7(e) are thus more robust than those obtained from Fig. 7(d) in case of errors of the velocity and direction measurements. If we reduce the noise level (e.g. *R* = 0.001), the ellipticity obtained with RTBF is tending to the true ellipticity even for the waves having the ellipticity with the smallest absolute value. No improvement can be observed with Poggi's method in the case of low noise by comparison with the results of Table 1.

We have shown in this section that the combination of the vertical component and the radial projection in the computation of the cross-spectrum matrix offers a real advantage over classical approaches used so far with single component or single projection processing for the determination of $\vec{k}$. We also provide a new way of extracting Rayleigh ellipticity in an efficient way. For a single far-field plane wave, the estimation is accurate even in the presence of strong incoherent noise. For multiple waves, the separation of 12 Rayleigh waves with only 12 three-component sensors in the presence of six Love waves is working fairly well to retrieve both the velocity and directions of propagation. For a 12 per cent portion of incoherent noise, about 5 per cent error is found on the ellipticity determination for waves with the strongest ellipticity (2). For the other waves with the lowest ellipticity (−0.5), the algorithm is not so efficient and larger errors are observed. The errors are directly linked to the level of incoherent noise. We also provide a way to quantify the level of incoherent noise that can be used to assess the reliability of the computed ellipticity values. In all cases, there is no ambiguity on the sign of the ellipticity. The CPU cost of the proposed method is still reasonable. For instance, 1 hr of ambient vibrations recorded on the vertical component is processed in 49 s on a 16-core machine. The three-component processing in the same conditions lasts approximately 380 s, about 7 to 8 times longer. In the following, we apply the new methodology on a synthetic and a real case.

## 3 SYNTHETIC AMBIENT VIBRATIONS

We test RTBF on model M2.1 commonly used for comparison of methods within SESAME European Project ('Site EffectS assessment using AMbient Excitations', Bard et al. 2004, Bonnefoy-Claudet et al. 2006) and in particular by Fäh et al. (2008), Poggi & Fäh (2010), Maranò et al. (2012) and Maranò et al. (2017) for three-component processing. The propagation medium is made of a single soft sediment layer on top of a rigid half-space, with a strong impedance contrast. The original data set provides 405 s of simulated noise up to 15 Hz over 38 three-component sensors with source to receiver distances varying from 150 to 700 m. Green's functions were computed with the code provided by Hisada (1994) and Hisada (1995) and they were convolved with delta-like impulse sources having random orientations and amplitudes. The complete description of the model can be found in the aforementioned references. We consider only two sensor geometries: *all* which contain all 38 receivers of the data set (200 m aperture) and *array 2* (14 sensors, 80 m aperture) which is used by Poggi & Fäh (2010) and Maranò et al. (2017). The geometries of these arrays can be found in Bard et al. (2004), page 9 and Poggi & Fäh (2010), Fig. 7, respectively.

To follow the recommendations of Capon (1969) about the number of non-overlapping blocks (*M*) for the computation of averaged cross-spectral density matrices (*M* > *N*), much longer signals are required to obtain a reliable dispersion curve around the frequency
of the resonance peak observed for model M2.1 (2 Hz). For instance, with blocks of 100 periods per block (50 s at 2 Hz) and 38 stations, the minimum duration is 1900 s at 2 Hz without block overlapping and $M=N$. Hence the same set of receivers and sources is simulated here with a longer duration, up to 6 hr and sampled at 200 Hz. The source spectrum is extended up to 50 Hz. The source distribution in time is simply revised by defining a probability of occurrence of 1% at every time period. On average, each source is fired around 2140 times during 6 hr.

Sets of 4N blocks are used to compute averaged cross-correlation matrices ($2N \times 2N$), with N being the number of three-component sensors. Smaller sets were tested and provide less accuracy in the determination of dispersion and ellipticity curves. The length of blocks is fixed to 100 periods, hence blocks are shorter in time at high frequency than those at low frequency. For a given frequency, if enough signal duration is available, the operation is repeated severable times to generate a maximum of 50 estimations of the wavefield parameters (direction, slowness, ellipticity angle, incoherent noise ratio and beam power), equally distributed along the total duration. We forbid any overlap between blocks but block sets can overlap. Neither spectral smoothing, amplitude normalization, frequency average nor noise decomposition is performed. For each frequency and each block set, all local maxima of the beam power are automatically identified using the procedure described in Section 2.6. Histograms are then built for each frequency, for dispersion and ellipticity curves.

In Fig. 9, we show the raw results calculated with VBF and RTBF on array *all* without any selection. While the fundamental mode dispersion curve is correctly reconstructed with a good accuracy from 2.3 to 20 Hz with both methods, the first six higher-order modes are better defined and separated by RTBF. They can be picked easily without any subjective expertise. Below 2.3 Hz, all methods are losing accuracy with a strong bias around 2 Hz, the resonance frequency of model M2.1, due to the vanishing intensity on the vertical component. A similar phenomenon is observed at 5 Hz on the first higher-order mode, in which ellipticity angle is changing from $-\pi/2$ to $\pi/2$, that is, also a vanishing of the Rayleigh wave intensity in the vertical component. A small bias is also visible on the RTBF curve for the first higher-order mode. The second higher-order mode apparently suffers from a similar lack of energy at around 6 Hz.

The first four modes in Fig. 9(b) are manually and roughly picked. All samples in which slowness falls within a 20 per cent range around the picked curves are selected. Thanks to the cleanliness of Fig. 9(b) there are almost no ambiguity during the selection process. The four selected modal dispersion curves are plotted in Figs 10(a), (d), (g) and (j). Their corresponding ellipticities and estimated incoherent noise ratios $\hat{R}$ (eq. 40) are also shown in Fig. 10. For the fundamental mode (Figs 10a–c), the estimated ellipticity slightly deviates from the true ellipticity between 3.6 and 5.2 Hz which is associated with an estimated noise ratio of above 3. The noise ratio also increases at high frequency with a related increase of the ellipticity uncertainty visible above 15 Hz. The first higher-order mode (Figs 10d–f) deviates from the theoretical curves at 3.7, 5 and 8.8 Hz and in parallel, the incoherent noise ratio rises above 3 for these frequencies. For the second and third higher-order modes the valid range for ellipticity is rather restricted. It matches relatively well the frequency range where values of $\hat{R}$ are below 3. Therefore a level of 3 for the estimated incoherent noise ratio appears to be a suitable threshold to select the good ellipticity values. We already observe in Section 2.6 that the level of incoherent noise influences the ellipticity measure especially for the case where several

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dispersion curves obtained by processing 6-hr long signals simulated for array *all* made of 38 receivers. (a) Dispersion curves retrieved with a single-component VBF. (b) Dispersion curves produced by RTBF. All local maxima of the beam powers are reported without any selection. The colour scale is normalized to the maximum value of each plot keeping the same ratio ($10^4$) between the maximum and the minimum of the scale. The red straight lines are the theoretical array limits as specified by Wathelet et al. (2008), from left to right: $k_{min}/2$ (plain line), $k_{min}$ (dashed line), $k_{max}/2$ (dashed line), $k_{max}$ (plain line). The dotted curves are the theoretical dispersion curves for the fundamental and the first higher-order modes.</figcaption></figure>

waves with significantly different ellipticities values are propagating together. The most affected are the ellipticities with the lowest absolute values, for example, the fundamental mode between 3.6 and 5.2 Hz in which absolute value is overestimated. The fundamental mode ellipticity is particularly well retrieved here on the flat part between 5.2 and 25 Hz because the high number of sensors allows a proper separation of the majority of the passing waves, reducing to a small quantity the uncaught waves considered as uncorrelated noise.

In Fig. 11 the dispersion and ellipticity curves are shown after a selection of all samples with $\hat{R}$ below 3. Outliers like those found for the fundamental mode above 10 Hz can be rejected by analysing each mode individually and recalling that one single ellipticity is allowed per frequency and per mode. Dispersion curves are retrieved accurately except at low frequency. This common feature can be observed on all modes drifting from the true curves when reaching their lowest frequency point, identified at 2.2, 3.9 and 7 Hz for the fundamental and the first two higher-order modes, respectively. In terms of wavenumber they correspond to $k_{min}$ divided by 3.9, 2.7 and 1.5 for the fundamental and the first two higher-order modes,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Estimated ellipticities and incoherent noise ratios $\hat{R}$ for the first four modes selected from the RTBF results in Fig. 9(b). (a–c) Dispersion, ellipticity curve and estimated incoherent noise ratio $\hat{R}$ for the fundamental mode. (d–f) First higher-order mode. (g–i) Second higher-order mode. (j–l) Third higher-order mode. Colour scales, theoretical curves and array limits are defined in the same way as in Fig. 9.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Dispersion (a) and ellipticity (b) curves obtained after a selection of all samples with an estimated incoherent noise ratio below 3. Colour scales, theoretical curves and array limits are defined in the same way as in Fig. 9.</figcaption></figure>

respectively. For all of them, either the vertical or the radial component is vanishing at their low frequency end point. The proposed criterion based on $\hat{R}$ is not sufficient to reject these biased points. Looking towards low frequencies, we can note for all modes that whenever $\hat{R}$ starts to increase, the dispersion and ellipticity curves start to drift from the true values. The inflection points on $\hat{R}$ curves are located at 2.3, 4.1 and 7 Hz for the fundamental and the first two higher-order modes, respectively. Hence for this case, considering the low frequency inflection point of $\hat{R}$ curves appears as a conservative option to select accurate dispersion and ellipticity curves below the classical array limits based on $k_{min}$.

The dispersion curves computed for array 2 with VBF and RBF (VRBF) by Poggi & Fäh (2010) are compared to the dispersion curves obtained with RTBF, using the original SESAME data set lasting only 6 min and 45 s. We run RTBF with the same parameters as for the 6 hr of signal described here above, except that a block overlap of 75 per cent is tolerated to circumvent the lack of data. A suboptimal efficiency is expected for high-resolution techniques but this is the price of a fair comparison with other methods available in the literature. All samples with $\hat{R}$ higher than 3 are filtered out. Outliers are not removed. In Fig. 12 we show RTBF results with the same display parameters as those used in Poggi & Fäh (2010, fig. 8c1). The low frequency inflection point of $\hat{R}$ curve is around 2.4 Hz. For the first high-order mode the inflection point is around 4.3 Hz (not shown). All slowness values above these cut-off frequencies are accurate for the fundamental mode and only slightly

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Results obtained with RTBF on array 2 from the original SESAME data set lasting 6 min and 45 s. (a) Dispersion and (b) $\hat{R}$ curves. The axis scales are the same as in Poggi & Fäh (2010, fig. 8c1). A linear colour map is also chosen as in Poggi & Fäh (2010) except that the maximum of the colour scales is set to one-tenth of the maximum density. Theoretical curves and array limits are defined in the same way as in Fig. 9.</figcaption></figure>

biased for the first higher-order mode. Slowness uncertainties for the fundamental mode are smaller with RTBF than with VRBF. High probability density are achieved with RTBF between 2 and 4 Hz, preventing the usage of a truly normalized linear colour scale like in Poggi & Fäh (2010). The aliasing features visible above 5 Hz with VRBF disappear with RTBF. The second and the third higher-order modes are not properly separated by VRBF while RTBF detects the second higher-order mode without ambiguity with the third higher-order mode.

With a wavefield decomposition (WaveDec), Maranò et al. (2017) computed the signed ellipticity for array 2, using also the original SESAME data set. Dispersion and ellipticity curves obtained with RTBF are shown in Figs 13(a) and (b) with the same display parameters as in Maranò et al. (2017, Figs 4c–f). Dispersion curves are retrieved in a similar way for the two methods. However with RTBF, the range extends to higher frequencies, the separation of the fundamental and the first higher-order mode between 4 and 5 Hz is effective, the second higher-order mode is detected, the uncertainties are about twice smaller [e.g. standard deviations are explicitly calculated for the fundamental mode in Fig. 13(a)], and the variation of the mean values versus frequency is smoother.

WaveDec and RTBF nicely delineate the fundamental mode ellipticity from 2.3 to 10 Hz (Fig. 13b), except when the horizontal
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Results obtained with RTBF on *array 2* from the original SESAME data set lasting 6 min and 45 s. (a) Dispersion, (b) ellipticity and (c) $\hat{R}$ curves obtained with RTBF. The axis scales are the same as in Maranò *et al.* (2017, Fig. 4). In panel (a) yellow dots with black error bars are the statistics calculated from the background probability density (RTBF). Colour scales and theoretical curves are defined in the same way as in Fig. 9.</figcaption></figure>

component is vanishing around 3.8 Hz. RTBF does not provide any samples in this region while WaveDec provides reliable values. The first higher-order mode is also well retrieved by the two methods. The second higher-order mode is approximately recovered by RTBF. Mean ellipticity curves are extracted for all modes above the inflection point of $\hat{R}$ curves and shown in Fig. 14 on a classical log-log ellipticity plot. WaveDec better follows the sharpness of the theoretical ellipticity at the trough of the fundamental curve and at the peak of the first higher-order mode. WaveDec also better fit the flat part of the fundamental mode. However, RTBF accurately retrieves the right flank of the main peak of the fundamental mode for all calculated ellipticity values while WaveDec overestimates the ellipticity below 3 Hz. In Fig. 14 the results from RayDec (Hobiger *et al.* 2009) are also compared with RTBF mean curves. A good match

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Absolute value of ellipticity curves calculated for *array 2* with various methods. The axis scales are the same as in Maranò *et al.* (2017, Fig. 4b). The black dotted curves are the theoretical ellipticity curves for the first three modes. The coloured circles are the mean values calculated with RTBF results (red, green and blue for the fundamental and the first two higher-order modes, respectively). The grey curves (mean and standard deviation with plain and dashed lines, respectively) are obtained with RayDec (Hobiger *et al.* 2009).</figcaption></figure>

is found from 2.4 to 3.3 Hz. Below this range, RayDec curves are biased. Comparing Fig. 14 and Poggi & Fäh (2010, fig. 8c2), we can observe similar results except that RTBF provides the ellipticity at a lower frequency, 2.4 Hz instead of 3 Hz. The values found with VRBF slightly overestimates the ellipticity between 3 and 3.5 Hz while RTBF appears to be unbiased in that range.

For a simulated wavefield containing also Love and body waves, we have shown in this section that RTBF is able to provide accurate results for dispersion and ellipticity curves above the resonance frequency with a better resolution than VBF. A bias at low frequency is a common issue with any f-k technique. For conventional f-k there are some widely admitted limits (e.g. Wathelet *et al.* 2008). For high-resolution techniques no precise rules exist to limit the interpretation for large wavelengths. The two selection criteria based on $\hat{R}$ are efficient on the M2.1 case either for the full data set *all* or for the subset *array 2*. With further tests on other synthetic cases and experimental data sets, they may help to go beyond the strict limitations $k_{min}$ proposed by Wathelet *et al.* (2008).

## 4 EXPERIMENTAL CASE: MIRANDOLA, ITALY

During *InterPACIFIC* project (Garofalo *et al.* 2016a,b), three sites in western Europe were densely investigated with ambient vibration arrays and borehole techniques. Among them, Mirandola (Italy) located on soft sediments is chosen here to test the capabilities of RTBF. The location, the array geometries and the acquisition parameters are described in Garofalo *et al.* (2016b). We use the data sets called *PCI* to *PC4*, with apertures from 30 to 810 m.

In Figs 15(a) and (b), we show the dispersion and ellipticity curves obtained with RTBF with the same processing parameters as in the previous section, except for the automatic picking of f-k peaks and for the block overlap. We accept only the f-k peaks with an amplitude above 50 per cent of the maximum amplitude at each frequency and time. For the largest aperture, we are forced
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Results obtained with RTBF on Mirandola data sets PC1 to PC4. (a) Dispersion curves, (b) ellipticity curves, (c) $\hat{R}$ for mode 1 and (d) $\hat{R}$ for mode(s) 2. In panels (a) and (b), the brightest colours represent the probability densities of the samples selected with a maximum $\hat{R}$ of 3. The faded colours in the background are the probability densities of all samples without any limit on $\hat{R}$. The small circles with error bars are the statistics calculated on the fundamental mode curves provided by the 14 participants of InterPACIFIC project (Garofalo et al. 2016b; Fig. 3). The higher-order mode curves provided by two InterPACIFIC teams (G03 and G12) are also shown below the fundamental. The straight black lines in panel (a) are the theoretical array limits for the largest data set PC4 as specified by Wathelet et al. (2008): $k_{min}/2$ (plain line), $k_{min}$ (dotted line). The colour scale is normalized to the maximum value of each plot keeping the same ratio ($10^4$) between the maximum and the minimum of the scale. The grey curves in panel (b) (mean and standard deviation with plain and dashed lines, respectively) are obtained with RayDec (Hobiger et al. 2009).</figcaption></figure>

to introduce an overlap of 50 per cent between blocks to generate results down to 0.65 Hz with only 2 hr of recordings.

The selection of the maximum $\hat{R}$ has little influence over the quality of the dispersion curve determination. The retrieved dispersion curve for the fundamental mode (labelled mode 1) compares very well with the results from InterPACIFIC project. Only two teams over 14 provided the higher-order modes (G03 and G12). Between 3 and 10 Hz we observe discrepancies between RTBF results, that suggest the presence of several modes, and the previous results. For the ellipticity, a large part of the fundamental mode is rejected at low frequency. The selection also slightly increases the mean values of the higher-order mode(s) (labelled mode(s) 2 because the experimental curve may be a compound of several higher-order modes). In Fig. 15(c), at low frequency for the fundamental mode, we note the same variation of $\hat{R}$ as in the synthetic case M2.1: a relatively low noise section that quickly increases towards low frequencies. The low frequency inflection point can be located around 0.8 Hz. Dispersion and ellipticity curves should be disregarded below that frequency. The peak of the H/V curve is around 0.75 Hz. Above 0.8 Hz RTBF and Raydec results are in good agreement.

## 5 CONCLUSIONS

From the assumption that the ambient wavefield is dominated by surface waves, we have investigated the properties of a beamforming based on a compound cross-spectral matrix built with vertical and radial components. For a single plane wave, we have shown analytically that a high-resolution beamformer can resolve the following wavefield parameters: the velocity, the direction of propagation, the sign and the absolute value of the Rayleigh ellipticity. In that simple case, the ellipticity values can be determined even in the presence of strong incoherent noise. The conclusions have been extended to multiple incident plane waves, except that the quality of the results are no longer independent of incoherent noise intensity. We also provide a way to quantify the level of the incoherent noise for each estimation of the wavefield parameters.

The RTBF is applied to a synthetic ambient wavefield, model M2.1 from SESAME project. The results are compared with single-component beamforming on a dense array made of 38 sensors and with long recordings (6 hr). RTBF offers a real advantage with a clear identification of all excited modes, especially at high frequencies where classical methods usually fail due to aliasing. However at low frequency in the neighbourhood of the resonance frequency or when any of the components is vanishing a bias in slowness and ellipticity values cannot be avoided like any other f-k methods. We have suggested two criteria based on $\hat{R}$ to select accurate samples. They are found to be efficient for the M2.1 case. With further tests on other synthetic cases and experimental data sets, they may help to go beyond the strict limitations $k_{min}$ proposed by Wathelet et al. (2008).

We have also processed the Mirandola data set (InterPACIFIC project). The retrieved fundamental dispersion curve is compatible with previous results. Higher-order modes that roughly match previous results are also clearly detected. Note that further improvements of RTBF results could be reached with a proper diagonal loading as proposed by Li et al. (2003). More detailed analysis of the relationships between the experimental dispersion and ellipticity curves would be required, which will be the subject of future works related to applications of RTBF to experimental data sets.
ACKNOWLEDGEMENTS

This work is part of the development of AVIOS (Ambient Vibration Integrated Online Service) supported and funded by EPOS-IP. We gratefully thank Dr M. Hobiger and an anonymous reviewer for their critical reading of the original manuscript and their constructive remarks. All tools developed in this paper are distributed under the GNU Public License v3 and they are available for download at http://www.geopsy.org. ISTerre is part of Labex OSUG@2020 (ANR10 LABX56).

A tutorial about RTBF is available at http://www.geopsy.org/wik i/index.php/Geopsy-fk.

REFERENCES

Bard, P.-Y., Bonnefoy-Claudet, S., Cornou, C., Fäh, D., Kristek, J. & Moczo, P., 2004. Simulation of seismic noise: report on parameter studies, SESAME report D12.09, last downloaded from http://sesame.geopsy.org on March 2018, Tech. rep.

Boaga, J., Cassiani, G., Strobbia, C.L. & Vignoli, G., 2013. Mode misiden- tification in Rayleigh waves: ellipticity as a cause and a cure, Geophysics , 78 (4), EN17–EN28.

Bonnefoy-Claudet, S., Cornou, C., Bard, P.-Y., Cotton, F. & Kristek, J., 2006. H/V ratio: a tool for site effect evalutation. Results from 1D noise simulations, Geophys. J. Int. , 167 , 827–837.

Boore, D.M. & Toksöz, M.N., 1969. Rayleigh wave particle motion and crustal structure, Bull. seism. Soc. Am. , 59 (1), 331–346.

Capon, J., 1969. High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE , 57 (8), 1408–1418.

Dunkin, J.W., 1965. Computation of modal solutions in layered, elastic media at high frequencies, Bull. seism. Soc. Am. , 55 (2), 335–358.

Fäh, D., Stamm, G. & Havenith, H.-B., 2008. Analysis of three-component ambient vibration array measurements, Geophys. J. Int. , 172 (1), 199–213.

Garofalo, F. et al. , 2016a. Interpacific project: comparison of invasive and non-invasive methods for seismic site characterization. Part ii: inter- comparison between surface-wave and borehole methods, Soil Dyn. Earthq. Eng. , 82 , 241–254.

Garofalo, F. et al. , 2016b. Interpacific project: comparison of invasive and non-invasive methods for seismic site characterization. Part i: intra- comparison of surface wave methods, Soil Dyn. Earthq. Eng. , 82 , 222– 240.

Henderson, H.V. & Searle, S.R., 1981. On deriving the inverse of a sum of matrices, SIAM Rev. , 23 (1), 53–60.

Hisada, Y., 1994. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths, Bull. seism. Soc. Am. , 84 (5), 1456–1472.

Hisada, Y., 1995. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths (part 2), Bull. seism. Soc. Am. , 85 (4), 1080–1093.

Hobiger, M., Bard, P.-Y., Cornou, C. & Le Bihan, N., 2009. Single station determination of Rayleigh wave ellipticity by using the random decrement technique (RayDec), Geophys. Res. Lett. , 36 (14), 0–4.

Hobiger, M., Le Bihan, N., Cornou, C. & Bard, P.-Y., 2012. Multicomponent signal processing for Rayleigh wave ellipticity estimation: application to seismic hazard assessment, IEEE Signal Process. Mag. , 29 (3), 29–39.

Hobiger, M. et al. , 2013. Ground structure imaging by inversions of Rayleigh wave ellipticity: sensitivity analysis and application to European strong- motion sites, Geophys. J. Int. , 192 (1), 207–229.

Li, J., Stoica, P. & Wang, Z., 2003. On robust capon beamforming and diagonal loading, IEEE Trans. Signal Process. , 51 (7), 1702–1715.

Maranò, S., Reller, C., Loeliger, H.-A. & Fäh, D., 2012. Seismic waves es- timation and wavefield decomposition: application to ambient vibrations, Geophys. J. Int. , 191 (1), 175–188.

Maranò, S., Hobiger, M. & Fäh, D., 2017. Retrieval of Rayleigh wave elliptic- ity from ambient vibration recordings, Geophys. J. Int. , 209 (1), 334–352.

Poggi, V. & Fäh, D., 2010. Estimating Rayleigh wave particle motion from three-component array analysis of ambient vibrations, Geophys. J. Int. , 180 (1), 251–267.

Roueff, A., Roux, P. & Réfrégier, P., 2009. Wave separation in ambient seismic noise using intrinsic coherence and polarization filtering, Signal Process. , 89 (4), 410–421.

Scherbaum, F., Hinzen, K.-G. & Ohrnberger, M., 2003. Determination of shallow shear wave velocity profiles in the Cologne/Germany area using ambient vibrations, Geophys. J. Int. , 152 , 597–612.

Schmidt, R., 1986. Multiple emitter location and signal parameter estima- tion, IEEE Trans. Antennas Propag. , 34 (3), 276–280.

Tokimatsu, K., 1997. Geotechnical site characterization using surface waves, in 1st International Conference on Earthquake Geotechnical Engineering , vol. 3 , pp. 1333–1368, Balkema, Rotterdam.

Wagner, G.S., 1996. Resolving diversely polarized, superimposed signals in three-component seismic array data, Geophys. Res. Lett. , 23 (14), 1837– 1840.

Wathelet, M., Jongmans, D., Ohrnberger, M. & Bonnefoy-Claudet, S., 2008. Array performances for ambient vibrations on a shallow structure and consequences over V s inversion, J. Seismol. , 12 (1), 1–19.

Woods, J.W. & Lintz, P.L., 1973. Plane Waves at Small Arrays, Geophysics , 38 (6), 1023–1041.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

3c-algebra-v1.pdf

Please note: Oxford University Press is not responsible for the con- tent or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be di- rected to the corresponding author for the article.

APPENDIX A: RADIAL CROSS-SPECTRAL DENSITY MATRIX FOR RAYLEIGH AND LOVE WAVES

We develop the analytical expression of the radial cross-spectral density matrix for a wavefield made of one Rayleigh and one Love plane waves. Replacing X$_{Rz}$ by ̂ X$_{Rh}$ in eq. (4) and from the definition of ̂ X$_{Rh}$ in eq. (11), we have

$$
\frac{1}{M} \sum_{m=1}^{M} \hat{X}_{R_{h},m} \hat{X}_{R_{h},m}^{*} = \frac{1}{M} \sum_{m=1}^{M} \left[ \alpha^2 X_{R_{h}} X_{R_{h}}^{*} + \alpha_{L}^2 X_{L} X_{L}^{*} \right. \\
\left. + \alpha_{E}^2 \eta_{E} \eta_{E}^{*} + \alpha_{N}^2 \eta_{N} \eta_{N}^{*} \right. \\
\left. + \alpha \alpha_{L} X_{R_{h}} X_{L}^{*} + \alpha \alpha_{L} X_{L} X_{R_{h}}^{*} + \dots \right] . \quad (A1)
$$

If M is sufficiently large, the Rayleigh and Love waves are not correlated with the random noise on horizontal sensors ( η$_{E}$ and η$_{N}$ ) and the unspecified terms in eq. (A1) can be neglected. Rayleigh and Love waves are also produced by random sources (i.e. random phases) that may be considered as uncorrelated. Hence the last two specified terms can be disregarded as well. We also assume that the average incoherent noise intensity is the same for all horizontal components and for all locations, that is,

$$
\frac{1}{M} \sum_{m=1}^{M} \eta_{E} \eta_{E}^{*} = \frac{1}{M} \sum_{m=1}^{M} \eta_{N} \eta_{N}^{*} = |N_{H}|^2 \mathbf{I}, \quad (A2)
$$

where | N$_{H}$ | 2 is the incoherent noise energy for any arbitrary hor- izontal direction θ . The coherent signals for this direction for the
Rayleigh and Love waves are, respectively,

$$ 
\begin{aligned}
|S_{Rh}|_{0}^{2} &= \alpha^{2} |S_{Rh}|^{2} \\
|S_{L}|_{0}^{2} &= \alpha_{L}^{2} |S_{L}|^{2}.
\end{aligned}
\quad (A3) $$

If the ratio of incoherent noise over the coherent signals for horizontal components is the same as for vertical components (R), we can write that

$$ |N_H|^2 = R(\alpha^2 |S_{Rh}|^2 + \alpha_L^2 |S_L|^2). \quad (A4) $$

From eq. (8), the cross-products become

$$ 
\begin{aligned}
\frac{1}{M} \sum_{m=1}^{M} \hat{X}_{Rh,m} \hat{X}_{Rh,m}^{*} &= \alpha^2 |S_{Rh}|^2 \boldsymbol{q}_0 \boldsymbol{q}_0^* + \alpha_L^2 |S_L|^2 \boldsymbol{q}_1 \boldsymbol{q}_1^* + |N_H|^2 \boldsymbol{I} \\
&= \alpha^2 |S_{Rh}|^2 (\boldsymbol{q}_0 \boldsymbol{q}_0^* + R \boldsymbol{I}) \\
&\quad + \alpha_L^2 |S_L|^2 (\boldsymbol{q}_1 \boldsymbol{q}_1^* + R \boldsymbol{I}) \\
&= \alpha^2 e_0^2 (1 + \tau_0^2) |S_{Rz}|^2 (\boldsymbol{q}_0 \boldsymbol{q}_0^* + R \boldsymbol{I}) \\
&\quad + \alpha_L^2 |S_L|^2 (\boldsymbol{q}_1 \boldsymbol{q}_1^* + R \boldsymbol{I}).
\end{aligned}
\quad (A5) $$

## APPENDIX B: MAXIMA OF $P'_{Rs}$ VERSUS THE PRESUMED ELLIPTICITY

In this section we calculate the maxima of $P'_{Rs}$ by checking the roots of its derivative versus the presumed ellipticity $e$.

$$ P'_{Rs} = C \frac{e^2}{A^2}, \quad (B1) $$

where $C$ is a constant that does not depend upon $e$,

$$ 
\begin{aligned}
A(e) &= \left( \frac{e^2}{e_0^2} + \alpha^2(1 + \tau_0^2) \right) - \left( \left( \frac{e}{e_0} + \alpha \right)^2 + \alpha^2 \tau_0^2 \right) \gamma, \\
\gamma(\bar{k}) &= \frac{|B(\bar{k} - \bar{k}_0)|^2}{2 + \frac{R}{N}}.
\end{aligned}
\quad (B2) $$

$$ \frac{d P'_{Rs}}{de} = C \frac{2eA^2 - 2e^2 A \frac{dA}{de}}{A^4}, \quad (B3) $$

where

$$ \frac{dA}{de} = \frac{2e}{e_0^2} - \frac{2}{e_0} \left( \frac{e}{e_0} + \alpha \right) \gamma. \quad (B4) $$

We first show that $A$ is never null and has a quadratic form.

$$ A = \frac{e^2}{e_0^2}(1 - \gamma) - 2 \frac{e}{e_0}\alpha\gamma + \alpha^2(1 + \tau_0^2)(1 - \gamma). \quad (B5) $$

The discriminant $\Delta$ is

$$ 
\begin{aligned}
\Delta &= 4 \frac{\alpha^2}{e_0^2} \gamma^2 - 4 \alpha^2 (1 + \tau_0^2) \frac{1}{e_0^2} (1 - \gamma)^2 \\
&= 4 \frac{\alpha^2}{e_0^2} \left[ (1 + \tau_0^2)(2\gamma - 1) - \tau_0^2 \gamma^2 \right].
\end{aligned}
\quad (B6) $$

Because $|B|^2$ is varying from 0 to 1, $\gamma$ is always less than 0.5. Hence the factor inside the bracket is always negative and $A$ has no root for $e \in \mathbb{R}$.

The derivative can be simplified to

$$ \frac{d P'_{Rs}}{de} = \frac{2Ce}{A^3}(\gamma - 1)\left[\frac{e^2}{e_0^2} - \alpha^2(1 + \tau_0^2)\right], \quad (B7) $$

which has roots at $e = 0$ and at $e = \pm|e_0||\alpha|\sqrt{1 + \tau_0^2}$. For a null $e$, $P'_{Rs}$ is also null while it is positive for any other value. Hence this root corresponds to a minimum. The other two roots correspond to maxima of $P'_{Rs}$ because the continuous function $P'_{Rs}$ tends to 0 when $e$ tends to $\pm\infty$.