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

HAL
open science

# Single-point dispersion measurement of surface waves combining translation, rotation and strain in weakly anisotropic media: theory

Le Tang, Heiner Igel, Jean-Paul Montagner

► To cite this version:

Le Tang, Heiner Igel, Jean-Paul Montagner. Single-point dispersion measurement of surface waves combining translation, rotation and strain in weakly anisotropic media: theory. Geophysical Journal International, 2023, 235, pp.24-47. 10.1093/gji/ggad199. insu-04155672

HAL Id: insu-04155672

https://insu.hal.science/insu-04155672

Submitted on 7 Jul 2023

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_2.png" /></figure>

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

# Single-point dispersion measurement of surface waves combining translation, rotation and strain in weakly anisotropic media: theory

Le Tang<sup>1</sup>, Heiner Igel<sup>1</sup> and Jean-Paul Montagner<sup>2</sup>

Accepted 2023 May 11. Received 2023 May 2; in original form 2022 November 28

## SUMMARY

A new approach is proposed for measuring the dispersion relation of surface waves in weakly anisotropic media using a single station, which consists of translational displacement and rotation or strain. The azimuth-dependent surface wave dispersion curve can be directly obtained by using the amplitude ratio of translational displacement to rotation or strain. Using observations from earthquakes from a variety of azimuths allows us to characterize the anisotropy of subsurface media. The amplitude ratio gives the localized azimuth-dependent dispersion relations of the formation beneath a receiver without requiring knowledge of the source or structure along the path. The theory shows that in anisotropic media the coupled quasi-Rayleigh wave and quasi-Love wave will result in elliptically polarized rotational motions. In addition, rotational motion information allows the determination of backazimuth even in general anisotropic media. We carry out numerical experiments, investigate the effects of noise and degree of anisotropy and demonstrate the potential for field studies.

**Key words:** Rotational seismology; Seismic anisotropy; Theoretical seismology.

## 1 INTRODUCTION

The study of seismic anisotropy has turned into a field of its own, relating material properties, stress-aligned heterogeneities or pore space, with dynamic phenomena such as mantle flow or fluid flow directions in crustal rocks (Anderson 1965; Forsyth 1975; Crampin 1977; Anderson & Dziewonski 1982; Nataf et al. 1984; Montagner & Tanimoto 1991; Legendre et al. 2021). A complete understanding of seismic anisotropy aids in characterizing the subsurface structure (e.g. Gupta 1973; Legendre et al. 2021) and earthquake disasters (Teanby et al. 2004). Previous studies on seismic anisotropy focus on translational displacement, including P-wave azimuthal anisotropy (Hess 1964; Francis 1969), P-wave polarization (Schulte-Pelkum et al. 2001; Al-Lazki et al. 2004; Mutlu & Karabulut 2011), shear wave splitting (Vinnik et al. 1989; Silver 1996; Crampin & Chastin 2003), surface wave azimuthal anisotropy and polarization analyses (Hess 1964; Forsyth 1975; Tanimoto & Anderson 1985; Montagner & Nataf 1986; Maupin 1989; Montagner & Tanimoto 1991).

Azimuthal anisotropy based on surface wave dispersion can be studied by determining the dispersion curves of different azimuth angles using single station methods, interstation methods, or array-based methods. The single station method has been used principally to measure fundamental mode dispersion curves (Levshin et al. 1972; Forsyth 1975; Ekström et al. 1997). Some of the more sophisticated approaches require knowledge of prior information of the source function or approximate earth structure (e.g. Brown et al. 2022) limiting the domain of application. The interstation or array-based method (Forsyth et al. 2005; Pedersen 2006) makes use of phase differences between seismograms measured at different stations to estimate phase velocity. This has been successfully applied to seismic ambient noise (Yao et al. 2006, 2008, 2010; Stehly et al. 2009; Legendre et al. 2021) combined with cross-correlation analysis (Shapiro & Campillo 2004; Shapiro et al. 2005). The dispersion measurement of the traditional single-station method or multistation method is the average velocity between the source (earthquake source or ambient noise fictional source) and the station, even if the average interstation dispersion can be related to the local phase velocity via an integral over the interstation arc (Legendre et al. 2021). Heterogeneity will also greatly affect results using these methods because the initially extracted dispersion curve is the average dispersion information along the surface wave propagation path.

Single-point anisotropy analysis of the polarization motion of surface waves (Park & Yu 1993; Tanimoto 2004) seeks to retrieve azimuthal variations of amplitude or amplitude ratios and obtain constraints on the anisotropic structure directly under the receiver. However,

<footnote n="1"><sup>1</sup>Department of Earth and Environmental Sciences, Ludwig-Maximilians University Munich, 80333 Munich, Germany. E-mail: ltang@geophysik.uni-muenchen.de</footnote>

<footnote n="2"><sup>2</sup>Institut de Physique du Globe de Paris, Université de Paris, 75238 Paris, France</footnote>
the polarization anomaly of a given station depends not only on the structure below the station but also on the structure seen by the wave in a significant part of the propagation path, unless the structure is uniformly anisotropic over a wide range (Maupin 2004).

Because of the limitations (heterogeneity effect and complexity of polarization) of the above surface wave anisotropy analysis, this paper attempts to find a method, characterizing the azimuthal anisotropy of subsurface structure without the influence of heterogeneity along the path or source effects.

The rapid development in rotational sensor technology (Pancha et al. 2000; Schreiber et al. 2005; Bernauer et al. 2021; Zembaty et al. 2021; Igel et al. 2021) allows us to exploit the potential application of rotational motions. Under the assumption of an isotropic medium or transversely isotropic (VTI) medium with vertical symmetry axis, only an S-related wave (including conversions) can generate rotational motions (Tang & Fang 2021a). This allows us to uniquely identify and separate the wavefield using polarization characteristics (Sollberger et al. 2018; Schmelzbach et al. 2018; Wassermann et al. 2022). Because the rotational component of an S-related wave is perpendicular to the wave propagation direction and should have the same phase as the corresponding translational displacement, the propagation azimuth of an earthquake can also be easily estimated using a single station (Igel et al. 2007; Wassermann et al. 2016). This even works for the direction of moving seismic sources (Yuan et al. 2021).

Fichtner & Igel (2009) gave the verification that the sensitivity kernels of the amplitude ratio (translational displacement and rotation) of surface waves or S waves attain large absolute values only in the vicinity of the receiver, but not the source, and such measurements may be used for local shear-wave speed tomography. The estimation of local phase velocity measurements from amplitude ratios has been verified with real observations and a variety of different rotation sensors (Igel et al. 2007; Wassermann et al. 2016; Keil et al. 2021; Fang & Tang 2021). Tang & Fang (2023) gave theoretical formulas of cross-correlation between displacement and rotation, and theoretically verified that the local phase velocity can also be extracted from amplitude ratios derived from ambient seismic noise.

Recently, Noe et al. (2022) numerically exploited the characteristics of body waves in a general anisotropic medium, demonstrating the feasibility to invert local anisotropic elastic parameters from rotational motions and strain by estimating the body wave velocity from amplitude ratios. Consequently, we can expect to use the translational displacement and rotation or strain to retrieve locally anisotropic parameters of subsurface structure from surface waves. To demonstrate this novel processing approach is the key goal of this study.

The paper is organized into two main sections. First, we establish the theoretical basis of the azimuth-dependent dispersion measurements from a single seismic station in weakly anisotropic media. Secondly, we numerically investigate the capability of extracting the azimuth-dependent dispersion curve from a single station and illustrate the effect of coupling terms. Furthermore, we compare our approach to the more classic analysis using phase effects.

## 2 THEORY

In Cartesian Coordinates, the strain tensor and rotation can, respectively, be defined as (Aki & Richards 2002):

$$ \epsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) \qquad (1) $$

$$ \Omega_i = \frac{1}{2}(u_{k,j} - u_{j,k}), \qquad (2) $$

where $u_{i,j}$ represents the partial derivative of translational displacement of i component with respect to j component. i, j and k vary among x, y and z.

The relation between stress and strain can be given as:

$$ \sigma_{ij} = C_{ijkl}\epsilon_{kl}, \qquad (3) $$

where $C_{ijkl}$ (i, j, k and l vary among x, y and z) is the fourth-order elastic tensor and $\sigma_{ij}$ represents the stress tensor. Considering the free surface boundary condition, where we assume that the surface is horizontal, the traction is equal to zero which indicates $\sigma_{iz} = 0$. Therefore,

$$ \sigma_{iz} = C_{izkl}\epsilon_{kl} = 0. \qquad (4) $$

Under the assumption that the media is only weakly anisotropic, the translational displacement of the anisotropic media can be written as a combination of Rayleigh and Love wave eigenfunctions in the reference isotropic media (Tanimoto 2004), whose expressions are given below:

$$ \mathbf{u} = a_L \mathbf{u}_L + a_R \mathbf{u}_R, \qquad (5) $$

where $\mathbf{u}_L$ and $\mathbf{u}_R$ are the Love and Rayleigh wave translational displacement of the reference isotropic medium, respectively, $a_L$ and $a_R$ are the coefficients to be determined from the stationarity of the Lagrangian. The dispersion relation and polarization of translational displacement in anisotropic media derived from the stationarity of the Lagrangian under the assumption of eq. (5) have been verified by comparison with other direct integral results (Tanimoto 2004). Therefore, we will derive the polarization of rotation and strain under this assumption.

In isotropic media, $\mathbf{u}_L$ and $\mathbf{u}_R$ have the form:

$$ \left\{ \begin{array}{l} \mathbf{u}_L = [-\sin \psi W(z), \cos \psi W(z), 0] e^{i\kappa_0(\cos \psi x + \sin \psi y) - i\omega t} \\ \mathbf{u}_R = [\cos \psi V(z), \sin \psi V(z), iU(z)] e^{i\kappa_0(\cos \psi x + \sin \psi y) - i\omega t}, \end{array} \right. \qquad (6) $$
where $\sin\psi^2 + \cos\psi^2 = 1$, $\psi$ is the azimuth of the direction of propagation measured from the X-axis. $\kappa$ is the modulus of the horizontal wavenumber vector $\vec{\kappa}$. $\kappa_{l0}$ and $\kappa_{r0}$ represent the wavenumber of Love waves and Rayleigh waves in isotropic media, respectively. $W(z)$ is the depth-dependent eigenfunction of the Love wave, $V(z)$ and $U(z)$ of the Rayleigh wave. $\omega$ is the circular frequency.

Substituting eqs (5) and (6) into the Lagrangian (Tanimoto 2004), we obtain:

$$ \mathcal{L}(u_i^*, u_i) = \omega^2 \int_{-\infty}^{0} \rho u_i^* u_i dz - \int_{-\infty}^{0} C_{ijkl} \varepsilon_{ij}^* \varepsilon_{kl} dz. \qquad (7) $$

By using the relations $\partial\mathcal{L}/\partial a_L = 0$, $\partial\mathcal{L}/\partial a_R = 0$ (Tanimoto 2004), the coefficients in the limit of weak anisotropy under $|A - B| \ge E$ can be expressed as:

For quasi-Love (qL) waves:

$$ (a_L, a_R) = \left( 1, \frac{E}{A - B} \right). \qquad (8) $$

For quasi-Rayleigh (qR) waves:

$$ (a_L, a_R) = \left( -\frac{E}{A - B}, 1 \right). \qquad (9) $$

A and B are, respectively, the qL and Rayleigh waves' eigenfrequency square and E is the coupling term. The detailed expressions of E, A and B can be found in the appendix of Tanimoto (2004). The translational displacement of qL wave can be expressed from eqs (5), (6) and (8) as:

$$ \mathbf{u}_{qL} = [-\sin \psi W(z) + \cos \psi T V(z), \cos \psi W(z) + \sin \psi T V(z), i T U(z)] e^{i \kappa_l (\cos \psi x + \sin \psi y) - i \omega t}. \qquad (10) $$

The displacement of qR waves can be expressed from eqs (5), (6) and (9) as:

$$ \mathbf{u}_{qR} = [\sin \psi T W(z) + \cos \psi V(z), -\cos \psi T W(z) + \sin \psi V(z), i U(z)] e^{i \kappa_r (\cos \psi x + \sin \psi y) - i \omega t}, \qquad (11) $$

where $T = \frac{E}{A-B}$. $\kappa_l$ and $\kappa_r$ represent the wavenumber of qL waves and qR waves in anisotropic media, respectively.

Eqs (10) and (11) are obtained by Tanimoto (2004), indicating that because of the effect of coupling, the polarization of Love waves in anisotropic media becomes elliptical instead of the linear particle motion in isotropic media. The particle motion of Rayleigh waves has a transverse component, which is no longer confined in the radial–vertical plane. To derive the dispersion relations between displacement and rotation or strain, in the next section, we first derive the polarization of rotation and strain combined with the free surface boundary condition and then obtain the dispersion formulas.

## 2.1 Isotropic media

### 2.1.1 Polarization and dispersion relation of rotation and strain in isotropic media

Under the assumption of isotropic layered media, there is no azimuth-dependent velocity variation. Combining the boundary conditions eq. (4) and displacement function eq. (6), we can derive the following equations of eigenfunctions:

$$ \left\{ \begin{array}{l} V'(0) = \kappa_{r0} U(0) \\ (2\mu + \lambda)U'(0) + \lambda\kappa_{r0}V(0) = 0 \\ W'(0) = 0, \end{array} \right. \qquad (12) $$

where $\mu$ and $\lambda$ are the Lamé parameters for isotropic media. $U'$, $V'$ and $W'$ represent the depth derivatives ($\partial_z|_{z=0}$) of the first-order eigenfunctions. Substituting eqs (12) and (6) into eqs (1) and (2), we obtain the polarization of rotation and strain at the surface ($z = 0$) in Cartesian coordinate (see Appendix C for detailed expressions). To better understand the polarization of rotation and strain, we project the polarization of eqs (6) and (C1)–(C4) onto the R–T coordinate system ($\cos\psi = 1$, $\sin\psi = 0$) of the propagation direction (R direction) shown in Fig. 1 with two identical vertical axes.

In this paper, we analyse the polarization of rotational components $\Omega_i(i = r, t, z)$ and the radial strain component $\varepsilon_{rr}$, which correspond to rotational seismometers and strain sensors or distributed acoustic sensing (DAS), respectively. After the coordinate projection [see section 5.4 of Chou & Pagano (1992) for the transformation relationship], the polarization of translational displacement, rotation and strain at the surface ($z = 0$) in the R–T coordinate system can be expressed as For Love waves:

$$ \left\{ \begin{array}{l} u_r = 0 \\ u_t = W \\ u_z = 0 \end{array} \right. \qquad (13) $$

$$ \left\{ \begin{array}{l} \Omega_r = 0 \\ \Omega_t = 0 \\ \Omega_z = i \kappa_{l0} W / 2 \\ \varepsilon_{rr} = 0. \end{array} \right. \qquad (14) $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Diagram of coordinate transformation.</figcaption></figure>

For Rayleigh waves:

$$ \left\{ \begin{array}{l} u_r = V \\ u_t = 0 \\ u_z = iU \end{array} \right. \qquad (15) $$

$$ \left\{ \begin{array}{l} \Omega_r = 0 \\ \Omega_t = \kappa_{r0} U \\ \Omega_z = 0 \\ \varepsilon_{rr} = i\kappa_{r0}V, \end{array} \right. \qquad (16) $$

where we omit the exponential term $e^{ikr - i\omega t}$. In the following derivation we also omit this term.

From a theoretical point of view, it is known that the translational displacement of the Rayleigh wave shows elliptical polarization in the R-Z plane and rotation shows linear polarization in the T direction. Both translational displacement and rotation of the Love wave show linear polarization in two orthogonal directions, respectively. Previous studies of surface waves with rotational components assume isotropic structure (Igel et al. 2007; Kurrle et al. 2010; Keil et al. 2021; Fang & Tang 2021; Tang & Fang 2021b, 2023), except for the analysis of body waves by Noe et al. (2022). In the following, we will demonstrate that in anisotropic media rotational motions are no longer linearly polarized.

Comparing eqs (13)–(16), the surface wave dispersion relations can be directly obtained by combining the polarization of translational displacement and rotation or strain.

For Love waves:

$$ c_L(\omega) = \left| \frac{\dot{u}_t}{2\Omega_z} \right|. \qquad (17) $$

For Rayleigh waves:

$$ \left\{ \begin{array}{l} c_R(\omega) = \left| \frac{\dot{u}_r}{\varepsilon_{rr}} \right| \\ c_R(\omega) = \left| \frac{\dot{u}_z}{\Omega_t} \right|, \end{array} \right. \qquad (18) $$

where $\dot{u}$ represents the first-order temporal derivative of displacement. Eqs (17) and (18) indicate that the dispersion phase velocity of Love and Rayleigh waves can directly be measured by the ratio of displacement to rotation or strain. This has been used widely in previous studies (Igel et al. 2007; Kurrle et al. 2010; Keil et al. 2021; Fang & Tang 2021).

## 2.2 Anisotropic media

Because of the coupling between qR and qL waves in (weakly) anisotropic media (Tanimoto 2004), the polarization of rotation of the two waves is no longer linear and shows elliptical form. An anisotropic medium is generally characterized by azimuth-dependent velocity variations. In the following derivation, we will separately consider three media: vertical transversely isotropic (VTI) media, general orthorhombic media [horizontal transversely isotropic (HTI) media is included] and general anisotropic media, whose expressions of the elastic tensor can be found in Appendix A. Among them, qL and qR waves are completely decoupled in VTI media and the phase velocity does not show azimuthal anisotropy. They are also decoupled for waves propagating along the symmetry axis and its perpendicular direction for HTI media.

### 2.2.1 Polarization and dispersion relation of rotation and strain in VTI media

The translational motion of VTI media has the same form as for isotropic media which can also be expressed by eq. (6). Therefore, the displacement, rotation, and strain also have the same polarization form, except for the vertical strain component $\varepsilon_{zz} = -\frac{c_{13}}{c_{33}}i\kappa V$ of qR wave, characterized by different elastic parameters. In addition, the dispersion formulas of VTI media generally also show the same form as isotropic media in eqs (17) and (18).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Diagram of polarization of translational displacement and rotation in isotropic media. The translational displacement of Rayleigh wave is elliptical polarization, while the love wave is linear polarization. The rotations of Love and Rayleigh waves are linearly polarized.</figcaption></figure>

### 2.2.2 Polarization and dispersion relation of rotation and strain in general orthorhombic media

Orthorhombic medium with nine independent parameters is characterized by azimuth-dependent velocity variations derived using the free surface boundary conditions eq. (4) as:

For qL waves:

$$ 
\left\{ 
\begin{aligned}
& \sin \psi T V'(0) + \cos \psi W'(0) = \kappa_l \sin \psi T U(0) \\
& - \sin \psi W'(0) + \cos \psi T V'(0) = \cos \psi \kappa_l T U(0) \\
& C_{33} T U'(0) + C_{13} \kappa_l \cos \psi [-\sin \psi W(0) + \cos \psi T V(0)] + C_{23} \kappa_l \sin \psi [\sin \psi T V(0) + \cos \psi W(0)] = 0.
\end{aligned}
\right. 
\qquad (19a) 
$$

For qR waves:

$$ 
\left\{ 
\begin{aligned}
& \sin \psi V'(0) - \cos \psi T W'(0) = \kappa_r \sin \psi U(0) \\
& \sin \psi T W'(0) + \cos \psi V'(0) = \cos \psi \kappa_r U(0) \\
& C_{33} U'(0) + C_{13} \kappa_r \cos \psi [\sin \psi T W(0) + \cos \psi V(0)] + C_{23} \kappa_r \sin \psi [\sin \psi V(0) - \cos \psi T W(0)] = 0.
\end{aligned}
\right. 
\qquad (19b) 
$$

We substitute eqs (19a), (19b), (10) and (11) into eqs (1) and (2), and the polarization of rotation and strain in Cartesian coordinates can, respectively, be found in Appendix C. We project the polarization of eqs (10), (11) and (C5)–(C8) onto the coordinate system (R–T) of propagation direction (R direction) shown in Fig. 1 and the polarization of translational displacement, rotation, and strain at the surface (z = 0) can be expressed:

For qL waves:

$$ 
\left\{ 
\begin{aligned}
& u_r = T V \\
& u_l = W \\
& u_z = i T U
\end{aligned}
\right. 
\qquad (20) 
$$

$$ 
\left\{ 
\begin{aligned}
& \Omega_r = 0 \\
& \Omega_l = \kappa_l T U \\
& \Omega_z = i \kappa_l W / 2 \\
& \varepsilon_{rr} = i \kappa_r T V.
\end{aligned}
\right. 
\qquad (21) 
$$

For qR waves:

$$ 
\left\{ 
\begin{aligned}
& u_r = V \\
& u_l = -T W \\
& u_z = i U
\end{aligned}
\right. 
\qquad (22) 
$$

$$ 
\left\{ 
\begin{aligned}
& \Omega_r = 0 \\
& \Omega_l = \kappa_r U \\
& \Omega_z = -i \kappa_r T W / 2 \\
& \varepsilon_{rr} = i \kappa_r V.
\end{aligned}
\right. 
\qquad (23) 
$$

As illustrated in Fig. 3, the translational polarization of the qR wave deviates from the vertical–radial plane in general orthorhombic media while the qL wave transforms into elliptical polarization, in which the coupled qR wave will appear on the transverse component and the coupled qL wave will appear on the vertical component (Tanimoto 2004). The rotational components of both waves exhibit elliptical polarization. However, both are only polarized in the vertical–transverse plane which is perpendicular to the direction of propagation.

Comparing eqs (20)–(23), the surface wave dispersion relations between translational displacement and rotation or strain can be easily obtained in general orthorhombic media:

For qL waves:

$$ c_L(\omega, \psi) = \left| \frac{\dot{u}_l}{2\Omega_z} \right| \qquad (24a) $$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Diagram of polarization of translational displacement and rotation in general orthotropic media. The translational displacement and rotations of Love and Rayleigh waves are elliptically polarized.</figcaption></figure>

$$ \begin{cases} c_L(\omega, \psi) = \left| \frac{\dot{u}_r}{\varepsilon_{rr}} \right| \\ c_L(\omega, \psi) = \left| \frac{\dot{u}_z}{\Omega_t} \right|. \end{cases} \qquad (24b) $$

For qR waves:

$$ \begin{cases} c_R(\omega, \psi) = \left| \frac{\dot{u}_r}{\varepsilon_{rr}} \right| \\ c_R(\omega, \psi) = \left| \frac{\dot{u}_z}{\Omega_t} \right|. \end{cases} \qquad (25a) $$

$$ c_R(\omega, \psi) = \left| \frac{\dot{u}_t}{2\Omega_z} \right|. \qquad (25b) $$

Eqs (24a) and (25a) give the local azimuth-dependent dispersion relation of general orthorhombic media. An interesting point is that the additional dispersion of eqs (24b) and (25b) for general orthorhombic medium come from the coupled waves contributions but they are generally smaller. From a theoretical point of view, the derived dispersion formulas from the coupling term contribution are correct. However, from an observational point of view, measuring the dispersion curves using this approach is prone to errors due to the amplitude of the coupling term. When the coupling term is equal to zero, $T = 0$, the eqs (24a), (24b), (25a) and (25b) degenerate into the isotropic dispersion equations.

### 2.2.3 Polarization and dispersion relation of rotation and strain in general anisotropic media

Under the assumption of first-order perturbations, general anisotropic media with 21 independent variables degenerate into 13 (monoclinic media in Appendix A; Smith & Dahlen 1973; Montagner & Nataf 1986). The free surface boundary conditions eq. (4) results in the following form:

For qL waves:

$$ \left\{ \begin{aligned} C_{45}[\cos \psi W'(0) + \sin \psi T V'(0) - \kappa_l \sin \psi T U(0)] + C_{55}[\cos \psi T V'(0) - \sin \psi W'(0) - \kappa_l \cos \psi T U(0)] &= 0 \\ C_{44}[\cos \psi W'(0) + \sin \psi T V'(0) - \kappa_l \sin \psi T U(0)] + C_{45}[\cos \psi T V'(0) - \sin \psi W'(0) - \kappa_l \cos \psi T U(0)] &= 0 \\ C_{33} T U'(0) = C_{13} \kappa_l \cos \psi [\sin \psi W(0) - \cos \psi T V(0)] - C_{23} \kappa_l \sin \psi [\cos \psi W(0) + \sin \psi T V(0)] - \frac{C_{36}}{2} [2 \sin \psi \cos \psi \kappa_l T V(0) + \cos \psi^2 \kappa_l W(0) - \sin \psi^2 \kappa_l W(0)]. \end{aligned} \right. \qquad (26a) $$

For qR waves:

$$ \left\{ \begin{aligned} C_{45}[-\cos \psi T W'(0) + \sin \psi V'(0) - \kappa_r \sin \psi U(0)] + C_{55}[\cos \psi V'(0) + \sin \psi T W'(0) - \kappa_r \cos \psi U(0)] &= 0 \\ C_{44}[-\cos \psi T W'(0) + \sin \psi V'(0) - \kappa_r \sin \psi U(0)] + C_{45}[\cos \psi V'(0) + \sin \psi T W'(0) - \kappa_r \cos \psi U(0)] &= 0 \\ C_{33} U'(0) = C_{23} \kappa_r \sin \psi [\cos \psi T W(0) - \sin \psi V(0)] - C_{13} \kappa_r \cos \psi [\cos \psi V(0) + \sin \psi T W(0)] - \frac{C_{36}}{2} [2 \sin \psi \cos \psi \kappa_r V(0) + (\sin \psi^2 - \cos \psi^2) \kappa_r T W(0)]. \end{aligned} \right. \qquad (26b) $$

Due to the complexity of the free surface boundary conditions which cannot simplify, the following derivation except for the third term of eqs (26a) and (26b), we directly substitute eqs (10) and (11) into eqs (1) and (2). The polarization in Cartesian coordinates can be found in Appendix C. We project the polarization of eqs (C9)–(C12) onto the coordinate system (R–T) of propagation direction (R direction) shown in Fig. 1, where translational components are the same as general orthotropic media, and the polarization of rotation and strain at the surface ($z = 0$) can be expressed as:

For qL waves:

$$ \left\{ \begin{aligned} \Omega_r &= 0 \\ \Omega_t &= (T V' + \kappa_l T U)/2 \\ \Omega_z &= i \kappa_l W/2 \\ \varepsilon_{rr} &= i \kappa_l T V. \end{aligned} \right. \qquad (27) $$
For qR waves:

$$ 
\left\{ 
\begin{array}{l}
\Omega_r = 0 \\
\Omega_t = (V' + \kappa_r U)/2 \\
\Omega_z = -i\kappa_r T W/2 \\
\epsilon_{rr} = i\kappa_r V.
\end{array}
\right. 
\qquad (28) 
$$

The rotational polarization form in general anisotropic media is the same as that of general orthotropic media shown in Fig. 3. The polarization plane is perpendicular to the propagation direction, which means that even in a generally anisotropic medium, we can still use the polarization information of the rotation to obtain the backazimuth. There is a difference with orthotropic media where it is easy to obtain the dispersion curves using $\dot{u}_z/\Omega_t$. The additional term $V'$ in $\Omega_t$ will affect the ratio value, this point is discussed later.

Comparing eqs (20), (22) with eqs (27) and (28), we can obtain the surface wave dispersion relation between displacement and rotation or strain:

For qL waves:

$$ c_L(\omega, \psi) = \left| \frac{\dot{u}_t}{2\Omega_z} \right| \qquad (29a) $$

$$ c_L(\omega, \psi) = \left| \frac{\dot{u}_r}{\epsilon_{rr}} \right|. \qquad (29b) $$

For qR waves:

$$ c_R(\omega, \psi) = \left| \frac{\dot{u}_r}{\epsilon_{rr}} \right| \qquad (30a) $$

$$ c_R(\omega, \psi) = \left| \frac{\dot{u}_t}{2\Omega_z} \right|. \qquad (30b) $$

The dispersion eqs (29b) and (30b) in general anisotropic media come from the coupled wavefield contributions that are generally smaller. Due to the complexity of the medium, the rotational polarization is difficult to simplify into a simple form in combination with the free surface boundary condition equation. Thus, there is no rotational formula that can be used to calculate the qR wave dispersion curve.

## 2.3 Effect of the coupling on the calculation of dispersion curves

The complete decoupling of Rayleigh and Love waves in isotropic and VTI media allows us to easily measure the local dispersion curves. On the contrary, the coupling effect in weakly anisotropic media results in polarization planes of the displacement not perpendicular or parallel to the direction of propagation. Furthermore, the translational motion of the qR wave is no longer limited to the vertical–radial plane and the qL wave is no longer limited to the radial–transverse plane. There is no clear separation between the two waves unless the group velocity difference is large. Consequently, we have to assess the effect of the coupling when we apply our theory to real data. The effect of the coupling term T depends on the phase velocity difference of the two waves and the degree of anisotropy. This can be described by explicit formulas (Tanimoto 2004):

$$ T = \frac{E}{A - B}. \qquad (31) $$

The numerator and denominator of eq. (31) are simultaneously divided by the square of the wavenumber $\kappa$:

$$ T = \frac{E/\kappa^2}{A/\kappa^2 - B/\kappa^2} = \frac{E/\kappa^2}{c_R^2 - c_L^2}. \qquad (32) $$

Eq. (32) indicates that if we fix the wavenumber, and the phase velocity difference of two waves is large enough, the T term tends to zero. Then the coupling term can be neglected. When we use eqs (24a) and (25a) to calculate the dispersion curves and we cannot identify the coupled qR and qL waves, using eqs (20)–(23), we obtain:

$$ \left| \frac{\dot{u}_z}{\Omega_t} \right| = \left| \frac{[\omega T U]_{qL} + [\omega U]_{qR}}{[\kappa_l T U]_{qL} + [\kappa_r U]_{qR}} \right| \qquad (33) $$

$$ \left| \frac{\dot{u}_t}{2\Omega_z} \right| = \left| \frac{[\omega W]_{qL} + [-\omega T W]_{qR}}{[\kappa_j W]_{qL} + [-\kappa_r T W]_{qR}} \right| \qquad (34) $$

As illustrated in eqs (33) and (34), the vertical translational displacement and transverse rotation contain not only the dominant qR wave energy, but also the qL wave signal generated by the coupling effect, and the transverse translation and vertical rotation contain not only the dominant qL wave energy but also the qR wave signal. As shown in eq. (32), the effect of the coupling term T will decrease and tends to zero when the velocity difference between the two modes increases. Consequently, eqs (33) and (34) degenerate to:

$$ \left| \frac{\dot{u}_z}{\Omega_t} \right| \approx \left| \frac{[\omega U]_{qR}}{[\kappa_r U]_{qR}} \right| = c_R \qquad (35) $$

$$ \left| \frac{\dot{u}_t}{2\Omega_z} \right| \approx \left| \frac{[\omega W]_{qL}}{[\kappa_j W]_{qL}} \right| = c_L. \qquad (36) $$
<table><caption>**Table 1.** Layer properties of Model 1-reference isotropic media.</caption><thead><tr><th>Layer</th><th>Depth (km)</th><th>V<sub>p</sub> (km s<sup>-1</sup>)</th><th>V<sub>s</sub> (km s<sup>-1</sup>)</th><th>ρ (kg m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>1</td><td>0–30</td><td>6.6</td><td>3.8</td><td>3000</td></tr><tr><td>2</td><td>30–∞</td><td>8.0</td><td>4.6</td><td>3300</td></tr></tbody></table>

It is generally appropriate to calculate the dispersion curves using eqs (35) and (36) for long-period teleseismic records because of the large velocity difference between fundamental qR wave mode and qL wave modes in that frequency range. The coupling strength is very weak regardless of whether the waveform can be distinguished in the time window. At this time, we do not need to separate the weakly coupled waveform. In contrast, when the velocity difference between qR and qL waves is small, the coupling effect will become stronger. At this time, it is difficult to distinguish the two coupled waves which will limit the validation of the derived dispersion formulas. In the numerical section, we will compute the effect of the coupling term on measuring dispersion curves for data processing, and demonstrate that if the radiation energy of the focal mechanism is large enough, it will weaken the effects of the coupling on the amplitude ratio.

## 2.4 Analytical azimuth-dependent dispersion relation

To verify the derived dispersion relations above, we use the following analytical solution, which is proposed by Smith & Dahlen (1973). Mochizuki (1986) and Tanimoto (1986) also showed the equivalent formulas in a spherical earth for the first time. We use these to benchmark our numerical results in the next section.

By solving eq. (7) simultaneously, the azimuth-dependent dispersion relation is obtained (Tanimoto 2004), having the same form as derived by Smith & Dahlen (1973); Montagner & Nataf (1986) after neglecting the coupling term E.

$$ \delta c_L(\omega, \psi) = \frac{1}{2c_{L0}(\omega)} [L_1(\omega) + L_2(\omega) \cos(2\psi) + L_3(\omega) \sin(2\psi) + L_4(\omega) \cos(4\psi) + L_5(\omega) \sin(4\psi)] \quad (37a) $$

$$ c_L(\omega, \psi) = \delta c_L(\omega, \psi) + c_{L0}(\omega) \quad (37b) $$

$$ \delta c_R(\omega, \psi) = \frac{1}{2c_{R0}(\omega)} [R_1(\omega) + R_2(\omega) \cos(2\psi) + R_3(\omega) \sin(2\psi) + R_4(\omega) \cos(4\psi) + R_5(\omega) \sin(4\psi)] \quad (38a) $$

$$ c_R(\omega, \psi) = \delta c_R(\omega, \psi) + c_{R0}(\omega). \quad (38b) $$

Here we define $\psi$ as the azimuth of the wavenumber vector measured anticlockwise from the X-axis. $\delta c_L(\omega, \psi)$ and $\delta c_R(\omega, \psi)$ are the first-order perturbations in phase velocity dispersion of azimuth frequency-dependent qL and qR waves. $c_{L0}(\omega)$ and $c_{R0}(\omega)$ are phase velocity of Rayleigh and Love waves, respectively, for a reference isotropic medium. $c_L(\omega)$ and $c_R(\omega)$ represent the phase velocity of the qL wave and qR wave, respectively. $L_i(\omega)$ and $R_i(\omega)$ ($i = 1, 2, 3, 4, 5$) are, respectively, depth integration functions that involve some elastic parameters and eigenfunctions, where we used a simple integration expression derived by Montagner & Nataf (1986), whose explicit expressions can be found in eqs (2), (4) and (5) of Montagner & Nataf (1986). We use the generalized reflection and transmission coefficients method (Chen 1993) to calculate the eigenfunctions which have been used successfully in many cases (Chen 1999; Tang & Fang 2021a).

It is known that eqs (37a), (37b), (38a) and (38b) are derived based on the assumption of first-order perturbation which can be called the ordinary perturbation method without the contribution of mode coupling. An implicit assumption when utilizing this solution is that the difference between any two-mode phase velocity is larger than their respective perturbations (Maupin 1989). This is generally valid for the fundamental mode. Consequently, we will focus on the analysis of the fundamental modes in the following numerical section.

# 3 NUMERICAL ANALYSIS

In this section, we test the theory presented above for different anisotropic models. We analyse the effect of surface wave coupling on the estimation of dispersion curves. Synthetic data are generated by simulating the complete seismic wavefield applying a 3-D standard staggered-grid finite difference method (Fang et al. 2014) which has been benchmarked with the generalized reflection and transmission coefficients method (Tang & Fang 2021a). The free surface boundary condition is implemented, and the perfectly-match-layer (PML) boundary condition is applied to the sides and bottom boundary. A moment tensor source with a Ricker wavelet source time function is used to generate the synthetic waveforms. We output translational acceleration components $A_i(i = r, t, z)$ and rotational velocity components $\Omega_i(i = r, t, z)$. For the sake of simplicity, the synthetic data generated in a two-layer model will be used to validate the applicability of the derived theoretical formulas for all models.

## 3.1 Reference isotropic case

Model 1 is set up to be the reference isotropic medium, whose parameters are listed in Table 1. This is used to calculate the eigenfunction W, V and U of eqs (6), (10) and (11) and in the analytical solution of eqs (37a) and (38a). We focus on long-period seismograms which means that
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Theoretical dispersion curves of the reference isotropic Model 1. Rayleigh-0th: The fundamental mode of Rayleigh wave. Rayleigh-1st: The first higher mode of Rayleigh wave. Love-0th: The fundamental mode of Love wave. Love-1st: The first higher mode of Love wave.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distributions of sources (black circles) and station (black triangle) of reference isotropic Model 1.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Beachball of focal mechanism in eqs (39) and (40).</figcaption></figure>

we will neglect the higher-mode surface waves. However, the numerical modelling has to ensure that there is no (or negligible) higher-mode signal in the selected period range. Then the selected period range is also used for the following weakly anisotropic models.

Fig. 4 shows that for periods longer than 12 s, there is only the fundamental (0th) mode solution for Model 1. Therefore, in the following numerical analysis, we focus on the period range of 12–80 s. As illustrated in Fig. 5, all sources (black circles) are located at (0, 0, 1 km), whose azimuth $\psi$ is evenly distributed at 15° ranging from 0° to 90°. A station (black triangle) is located at the surface. The radial distance between all sources and the station is equal to 2500 km. The source–receiver geometry of Fig. 5 is used for the following models to verify the effectiveness of measuring the azimuth-dependent dispersion characteristics in anisotropic media.

Two sets of synthetic data with different central frequencies, 0.015 Hz (period = 66.7 s) and 0.034 Hz (period = 29.4 s) are generated to ensure a 12–80 s wide-band seismogram. The magnitude of the moment tensor source is given as (Beachball is shown in Fig. 6a)

$$ M = \begin{bmatrix} 0.1 & 0.2 & 0.35 \\ 0.2 & 0.25 & 0.5 \\ 0.35 & 0.5 & 0.15 \end{bmatrix} (\text{N} \cdot \text{m}). \qquad (39) $$
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Normalized seismograms of different azimuths with a 0.034 Hz (period = 29.4 s) central frequency of reference isotropic Model 1. $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Normalized seismograms of different azimuths with a 0.015 Hz (period = 66.7 s) central frequency of reference isotropic Model 1. $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

In this model, the grid spacing is 5 and 3 km for the two different central frequencies, respectively, and the time increment is 50 ms for the 3-D finite difference simulations. Figs 7 and 8 show the seismograms using the same focal mechanism with different central frequencies, where translational acceleration $A_z$ and rotational velocity $\Omega_t$ are associated with Rayleigh waves and $A_t$ and $\Omega_z$ are associated with Love waves. This indicates that the radiation energy of the Rayleigh wave is evenly distributed in seven directions and the radiation energy of directions of S4–S6 of the Love wave is much smaller. In combination with eqs (17) and (18) and simulated seismograms of Figs 7 and 8, we can calculate the dispersion curves in different directions, where we introduce a robust weighted least-squares method based on time-frequency analysis to measure the ratio (see Appendix B). The weight function is equal to 1 for all data points in Model 1.

As illustrated in Fig. 9, the measured dispersion points of both the Rayleigh wave (red) and Love wave (blue) using the ratio dispersion of eqs (17) and (18) match with the analytical solution of the fundamental mode in any direction even with smaller radiation amplitude (see
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparisons of dispersion curves of different azimuths calculated using translation-rotation dispersion eqs (17) and (18) with seismograms in Figs 7 and 8, and analytical solution of the fundamental mode of reference isotropic Model 1. There is no correlation between the calculated results and the azimuth-dependent radiation pattern.</figcaption></figure>

<table><caption>**Table 2.** Layer properties of Model 2-HTI media.</caption><thead><tr><th>Layer</th><th>Depth (km)</th><th>ρ (kg m<sup>-3</sup>)</th><th>α<sub>V</sub> (km s<sup>-1</sup>)</th><th>α<sub>H</sub> - α<sub>V</sub> / α<sub>H</sub></th><th>β<sub>V</sub> (km s<sup>-1</sup>)</th><th>β<sub>H</sub> - β<sub>V</sub> / β<sub>H</sub></th><th>η</th></tr></thead><tbody><tr><td>1</td><td>0–30</td><td>3000</td><td>6.6</td><td>0</td><td>3.8</td><td>0</td><td>1</td></tr><tr><td>2</td><td>30–∞ (HTI)</td><td>3300</td><td>8.0</td><td>5 per cent</td><td>4.6</td><td>5 per cent</td><td>0.7721</td></tr></tbody></table>

source 4–5 of Love waves). In addition, it also indicates that there are no higher mode seismograms in the selected periods. Consequently, we will select the period range of 12–80 s to analyse the effect of coupling waves on dispersion measurements in various anisotropic models.

## 3.2 Contamination through coupled waves

Model 2 is set up to investigate the effect of coupling between qR and qL waves on calculating dispersion curves. For simplicity, we use a HTI medium with a symmetric axis parallel to the X-axis. Model 2 consists of a half-space HTI medium based on Model 1, where the first layer is the same as that of the reference isotropic Model 1. Table 2 lists the medium properties of Model 2. We use five independent parameters as defined by Takeuchi & Saito (1972) to describe the HTI medium of Model 2, where α<sub>V</sub> is the P-wave velocity along the symmetric axis, α<sub>H</sub> is the P-wave velocity perpendicular to the symmetric axis, β<sub>V</sub> is the S-wave velocity along the symmetric axis and β<sub>H</sub> is the SH-wave velocity perpendicular to the symmetric axis. The anisotropic strength of HTI medium for both P wave and S wave is also 5 per cent. The generated seismograms of Figs 10 and 11 of Model 2 are based on the same source–receiver geometry (Fig. 5), focal mechanism and period ranges of Model 1.

Combining dispersion eqs (24a) and (25a) of HTI media with seismograms in Figs 10 and 11, the dispersion curves are evaluated utilizing the least-square algorithm (Appendix B), where all data points are included indicating that the weight function is equal to 1. As illustrated in Fig. 12, the qR dispersion curves (red points) calculated by eq. (25a) generally match well with the theoretical dispersion curves (black points) showing azimuthal anisotropy. The phase velocity increases from 0 degrees to 90° in the long period range corresponding to HTI medium. The results (blue points) of qL wave in the period range of 12–30 s also follow the trend of theoretical solution due to the reason that the short-period wavefield propagates in the isotropic medium of the first layer without the effect of anisotropy. On the contrary, the results of the qL wave deviate from the theoretical solution in the period range of 35–80 s, especially in the direction of S4–S6, indicating that the qL wave is contaminated seriously by the coupled qR waves. As illustrated by the seismograms in Figs 10 and 11, the radiation amplitude of the qL wave in the direction of S4–S6 is much smaller, especially the S4 and S5 directions, so that the coupled qR waves cannot be ignored (see eq. 34). This has an important influence on the results. The radiation amplitude of qR waves in all directions is relatively large which can overshadow the influence of coupled qL waves.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Normalized seismograms of different azimuths with a 0.034 Hz (period = 29.4 s) central frequency of HTI Model 2 using the focal mechanism of eq. (39). $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Normalized seismograms of different azimuths with a 0.015 Hz (period = 66.7 s) central frequency of HTI Model 2 using the focal mechanism of eq. (39). $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

To further illustrate the fact that the coupled wavefield is non-negligible when the radiation amplitude is very weak, we only change the focal mechanism, leaving all other simulation conditions unchanged, using the moment tensor source can be expressed as (Beachball is shown in Fig. 6b)

$$M = \begin{bmatrix} 0.1 & 0.2 & -0.35 \\ 0.2 & -0.25 & 0.5 \\ -0.35 & 0.5 & 0.15 \end{bmatrix} (\text{N} \cdot \text{m}) \qquad (40)$$

As shown in Figs 13 and 14, the radiation energy of the qL and qR waves have strong amplitudes in all directions. We use the same dispersion equation as in Fig. 12 to calculate the dispersion curve, and the results are shown in Fig. 15. It can be seen that the dispersion curve of the qL wave has been improved compared with the results in Fig. 12. This is consistent with the theoretical dispersion trend, but there are still some deviations affected by the coupling. We further note that the deviation of the qL dispersion curve in Fig. 12 is mainly because the weak radiation amplitude will further amplify the coupling.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparisons of dispersion curves calculated using translation-rotation dispersion eqs (24a) and (25a) with seismograms in Figs 10 and 11, and analytical solution of the fundamental mode of HTI Model 2.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Normalized seismograms of different azimuths with a 0.034 Hz (period = 29.4 s) central frequency of HTI Model 2 using the focal mechanism of eq. (40). $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

### 3.3 Azimuth-dependent coupling effects

As the anisotropy of the medium increases, the influence of the coupled wavefield on the amplitude will be larger. However, it can be seen from Figs 13 and 14 that even though the radiation amplitude of S4–S6 of the qL wave is relatively large, the calculated dispersion curve shown in Fig. 15 still has a small deviation in certain azimuths. The difference in the strength of the azimuthal anisotropy of the two wave types leads to the coupling being a function of azimuth. Fig. 16 shows the *T* (coupling term) value related to the coupled wavefield calculated theoretically for wavenumber $\kappa = 5 \times 10^{-5}$ [the expression of *T* can be found in appendix of Tanimoto (2004)] in Model 2 which shows the correlation with angle. With this wavenumber, the coupling of S4–S5 is stronger, and the effect on the ratio value is greater. Fig. 17 shows the dispersion curve calculated under different degrees of anisotropy while other model conditions are kept constant. We conclude that the influence of the coupled wavefield in the direction of 30–75° (S3–S6) is more stronger.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Normalized seismograms of different azimuths with a 0.015 Hz (period = 66.7 s) central frequency of HTI Model 2 using the focal mechanism of eq. (40). $A_z$: translational acceleration in the vertical direction; $A_t$: translational acceleration in the transverse direction; $\Omega_t$: rotational velocity in the transverse direction; $\Omega_z$: rotational velocity in the vertical direction.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Comparisons of dispersion curves calculated using translation-rotation dispersion eqs (24a) and (25a) with seismograms in Figs 13 and 14, and analytical solution of the fundamental mode of HTI Model 2.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Azimuth-dependent coupling term $T$ value with wavenumber $\kappa = 5 \times 10^{-5}$ in HTI Model 2. The coupling strength is the largest at the azimuths of S4 and S5, and is decoupled in the directions of S1 and S7.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 17. Comparisons of dispersion curves calculated using different anisotropic parameters and analytical solutions of the fundamental mode of HTI Model 2.</figcaption></figure>

## 3.4 Seismogram stacking to stabilize solution

In order to weaken the effect of coupled wavefield on the results in some directions when the station is in the radiation node of the source, associated with a small amplitude, and to obtain a stable and reliable dispersion curve, we consider fitting data with multiple source seismograms in the same direction to obtain a reliable phase velocity. This is also in line with the consideration of actual data processing.

The media parameters and geometry of Model 3 are the same with that of model 2, except that 23 seismograms are generated in each direction with 23 different sources to estimate the dispersion value. The number of events is selected randomly, even one event with good radiation amplitude can lead to satisfactory results (see Fig. 15). In Model 3, the depth of sources in each direction is randomly distributed in the range of 5–15 km, the centre period of the wavelet is randomly distributed in the 30–80 s, and each magnitude of the focal mechanism is randomly selected from −1 to 1. Unlike in Models 1 and 2, we also calculate the radial strain seismogram $e_{rr}$ and radial translation seismogram $A_r$, which are associated with the qR wave, to illustrate the applicability of the translation-strain dispersion in eq. (25a).

The least-square method of time-frequency analysis in Appendix B is also performed to estimate the dispersion value. To reduce the influence of the coupled wavefield, the data with small radiation amplitude should be removed, which is demonstrated in the analysis of Model 2. A simple criterion is defined that the data points that are less than 10 per cent of the maximum amplitude in each periodic signal will be removed, and only more than 10 per cent of the energy is retained. This criterion is not fixed and can be adjusted according to the strength of actual anisotropy and noise level. Consequently, the weight function in Appendix B, which can also be called the filter function, can be expressed as:

$$wf(\omega, t_i) = \begin{cases} 0; & \left| \frac{A(\omega, t_i)}{\max[A(\omega, t_i)]} \right| < 10 \text{ per cent} \\ 1; & \text{otherwise.} \end{cases} \qquad (41)$$

Fig. 18 shows the time-period spectra of the first five seismograms (S1 direction) and the time-period spectra filtered by the weight eq. (41). It can be seen that after the processing of the weight function, the energy below 10 per cent in each period is directly removed. And the time-period spectra of seismograms associated with other directions (sources) can also be obtained in the same way, although we do not show these spectra. This filtering is also suitable for the suppression of random noise. We will discuss this effect of random noise in the next section. Fig. 19 shows the calculated dispersion results using the filtered spectra in Figs 18 combined with the dispersion eqs (24a) and (25a), indicating that the results are consistent with the theoretical solution, whether it is a qR wave (red points) or a qL wave (blue points). At the same time, the results (green points) in Fig. 19 obtained by the translation-strain dispersion eq. (25a) also match well the theoretical values. It demonstrates that by using multiple seismograms, the unfavourable influence of the coupled wavefield can be eliminated, resulting in a stable and reliable dispersion curve. These simulation results from Model 1–3 demonstrate the correctness and applicability of the derived dispersion formulas which also means that the ratio of translational displacement to rotation or strain can identify the azimuthal anisotropy from a single station 6C (three displacement components and three rotation components) measurement.

## 3.5 The effect of random noise

From the perspective of seismological observations, the influence of noise or the uncertainties of amplitude measurements cannot be ignored. In particular, this paper uses frequency-dependent amplitude information to obtain the dispersion curves, and small amplitudes perturbations may have a strong impact on the results (Kurrle et al. 2010). Therefore, the least-square solution based on linear regression in Appendix B to fit the data helps to get a stable solution. The multiple seismograms of qR waves in the direction of azimuth $\psi = 0$ (S1) in Model 3 are used as input data for our noise analysis. Fig. 20 shows the first five seismograms ($A_z$) of Model 3 perturbed with random Gaussian noise
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Period-normalized time-period spectra of the first five seismograms (S1 direction). From top to bottom: time-period spectra of $A_z$; filtered time-period spectra of $A_z$ using weight function of eq. (41); removed time-period spectra of $A_z$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Comparisons of dispersion curves calculated using translation-rotation and translation-strain dispersion eqs (24a) and (25a) and analytical solution of the fundamental mode of Model 3.</figcaption></figure>

with varying signal-to-noise ratio (SNR). It can be seen that when the SNR is equal to 2, the noise covers the whole time window, and the amplitude of each waveform will be greatly changed.

We add Gaussian random noise with different SNRs to $A_z$ and $\Omega_t$ at the same time. To simulate the influence of noise when processing multiple seismograms for one azimuth, each SNR is generated 100 times and we calculate the kernel density function (Botev *et al.* 2010) of the phase velocity distribution, which represents its probability density distribution. When calculating the phase velocity using the approach described in Appendix B, we still use the weight function in eq. (41) to exclude amplitudes less than 10 per cent in each period. Fig. 21 shows the kernel density function under different SNRs, and the black solid line is its theoretical dispersion curve. When the SNR is equal to 2, it can be seen that most of the energy is confined to the range within the 1 per cent error bar, especially between 40 and 60 s, benefiting from its large amplitude energy. For this period range 70–80 s with small amplitude energy, its density function energy distribution is more dispersed than that of other periods. With an improvement of the SNR, the energy of its kernel density function becomes closer to the theoretical
<figure><img src="image_23.png" /><figcaption>**Figure 20.** The first five seismograms ($A_z$, S1 direction) of Model 3 with different SNR.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Comparisons of dispersion probability energy distribution with different SNR, while each SNR is simulated for 100 times, and analytical solutions of the fundamental mode of HTI Model 3.</figcaption></figure>

dispersion curve. When the SNR is greater than 10, about 90 per cent of the energy is distributed between a 1 per cent error bar, showing that the dispersion curve obtained from this SNR lets us resolve anisotropy of 3–4 per cent in the medium.

Consequently, for real data processing, the seismic data with a high SNR and large radiation amplitude should be selected as much as possible to derive the dispersion value, and the solution obtained by using the linear regression method is more reliable.

## 3.6 Estimating local dispersion characteristics

The commonly used method to study anisotropy based on phase difference measurement dispersion is to extract the average velocity between the source and the station (earthquake data), or the average velocity between stations (seismic ambient noise data). Therefore, if the measured dispersion curve shows angular anisotropy, it is difficult to judge whether it is caused by heterogeneity or anisotropy.

Model 4 is used to demonstrate the advantage of the ratio method on lateral resolution compared with the traditional phase difference method in analysing azimuthal anisotropy. As shown in Fig. 22, the medium in the outer circle where the source is located in an isotropic medium, with parameters shown in Table 1. The inner cylinder with a radius of 500 km below the station is an HTI medium with 3 per cent anisotropy of the first layer of body waves, and the second layer is also a half-space HTI medium model with 4 per cent anisotropy whose parameters are shown in Table 3. The seismograms for Model 4, which are associated with qR waves, are generated using the moment tensor source of eq. (40) whose central frequencies are 0.034 Hz (period = 29.4 s) and 0.015 Hz (period = 66.7 s), respectively. We still
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Distributions of sources and receiver of Model 4.</figcaption></figure>

<table><caption>**Table 3.** Layer properties of the inner cylinder of Model 4.</caption><thead><tr><th>Layer</th><th>Depth (km)</th><th>ρ (kg m<sup>-3</sup>)</th><th>α<sub>V</sub> (km s<sup>-1</sup>)</th><th>α<sub>H</sub> - α<sub>V</sub> / α<sub>H</sub></th><th>β<sub>V</sub> (km s<sup>-1</sup>)</th><th>β<sub>H</sub> - β<sub>V</sub> / β<sub>H</sub></th><th>η</th></tr></thead><tbody><tr><td>1</td><td>0-30 (HTI)</td><td>3000</td><td>6.6</td><td>3 per cent</td><td>3.8</td><td>3 per cent</td><td>0.8489</td></tr><tr><td>2</td><td>30-∞ (HTI)</td><td>3300</td><td>8.0</td><td>4 per cent</td><td>4.6</td><td>4 per cent</td><td>0.8090</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Figure 23.** Comparisons of dispersion curves calculated using translation-rotation ratio eqs (24a) and (25a) and analytical solution of the fundamental mode of the inner cylinder in Model 4.</figcaption></figure>

use the least-square method in Appendix B to calculate the dispersion points, where the weight function is equal to 1 for all data points. As illustrated in Fig. 23, the estimated dispersion is consistent with the theoretical dispersion curve of the inner cylindrical medium, reflecting the surface wave dispersion characteristics of the medium directly below the station. The derived ratio dispersion equation gives the localized azimuthal-dependent dispersion relations of the formation right beneath a receiver. This suggests that the smaller the wavelength, the smaller the influence of heterogeneity. Therefore, we expect that the ratio dispersion method will provide a higher lateral resolution imaging result on azimuthal anisotropy analysis, where the lateral resolution can be smaller than one wavelength (Tang & Fang 2023).

## 3.7 Estimation of the azimuth using horizontal rotation

Accurate calculation of the azimuth of the wavefield incident on the station is crucial in the study of azimuthal anisotropy. It is the most basic parameter in the subsequent correct inversion of the anisotropy and the evaluation of the azimuth of mantle flow or fracture orientation. The error of the angle will directly lead to errors of the anisotropy parameters. As illustrated in Figs 2 and 3, the polarization characteristics of the
<figure><img src="image_27.png" /><figcaption>Figure 24. Comparisons of azimuth calculated using ratio in eq. (42) and analytical solutions with seismograms of Model 2 and Model 3.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25. Diagram of rotation angle $\psi_0$.</figcaption></figure>

rotation allow us to use the rotation component to calculate the incident azimuth of the phase velocity without distinguishing the properties of the medium and considering the effect of the coupled wavefield. The expression for calculating the azimuth can be given as

$$ \tan(\psi) = -\frac{\Omega_x}{\Omega_y}. \qquad (42) $$

The rotational seismograms of Model 2 and 3 with different focal mechanisms are used to verify the feasibility of using only the rotational component in eq. (42). Fig. 24 shows that the azimuth calculated from the rotational components is consistent with the theoretical azimuth and is not affected by the coupled wavefield despite the weak radiation amplitude S4–S6 in Figs 9 and 10 and contamination by coupling. The results, which are calculated with multiple seismograms in Model 3, are in good agreement with the theoretical values, indicating the applicability of eq. (42).

## 3.8 Anisotropy study based on the rotation angle

The polarization characteristics of rotations in Fig. 3 show that there is an additional angle between the propagation direction of the qR wave and the horizontal translation polarization direction in anisotropic media (see Fig. 25), which is called the rotation angle here. Similarly, there is also an additional angle between the quasi-Love wave and T-axis. This rotation angle is also a function of the propagation direction, so it can also be used to study azimuthal anisotropy. Taking the qR wave as an example, it can be seen from eq. (22) that the formula for calculating the rotation angle of the qR wave is

$$ \tan(\psi_0) = \frac{u_t}{u_r} = -\frac{W}{V}T, \qquad (43) $$

where the calculation of $u_t$ and $u_r$ from $u_x$ and $u_y$ requires the azimuth which can be calculated using rotation components in eq. 42. The tangent function of this angle is represented by the coupling term $T$, which is an azimuth-dependent function. Its calculation is similar to eq. (20) in Tanimoto (2004) who used the ratio of the vertical translation component of the coupled qL wave to the horizontal translation component of the qL wave to analyse anisotropy, except that the amplitude is different.

Eqs. (42) and (43) provide an alternative method for studying anisotropy. The advantage is that the azimuth estimation based on the rotation component will not be affected by the coupled wavefield, while the variables estimated based on the translation component will be affected by the coupled wavefield.
## 4 DISCUSSION

We derived concise polarization form of rotation and strain based on the first-order perturbation. This form allows us to clearly describe the particle motion and the dominant surface waves can be verified by comparing with the numerical results. However, it is very hard to verify the dispersion relation (eqs 24b and 25b) of the coupled waves whose amplitude is generally pretty small in a weakly anisotropic medium.

This first-order perturbation theory shows that the rotational vector is orthogonal to the wavenumber vector, so the radial rotational component $\Omega_r$ is always equal to zero. From the observational point of view, the recorded seismogram will contain scattering waves from all directions due to the heterogeneity in the Earth, especially in the short period range. Therefore, the signal in the radial rotational component of an earthquake comes from the scattering waves rather than from the effects of anisotropy. As the frequency decreases, the earth structure generally will become more homogeneous, the amplitude of the radial rotational component will tend to be zero. This amplitude-dependent dispersion measurement approach requires very high accuracy for amplitude measurements. Even a very small amplitude perturbation can results in larger deviation (Kurrle et al. 2010). Consequently, it is very necessary to use multiple seismogram to obtain a reliable velocity.

For a general anisotropic medium, as illustrated by eqs (30a) and (30b), there is no available rotational component that can be used to calculate the Rayleigh wave velocity because the retrieval of phase velocity from eq. (30b) which is from coupled waves is extremely difficult from real data. Considering the fact that the earth structure is a general anisotropic medium, and that the amplitude recorded by a rotational seismometer is more accurate than those of DAS at a single point observation, and the amplitude of fundamental Rayleigh waves is less affected by coupled waves than Love waves. So, is it possible to combine the rotational component to calculate the velocity of qR wave for a general anisotropic medium? Let us consider the assumption that $\kappa_{r0}$ is close to $\kappa_r$, it means that the degree of anisotropy is very weak. Then, $V' \approx \kappa_r U$ in eq. (28) because in isotropic media, $V' = \kappa_{r0}U$ (see eq. 12). Combining eqs (28) and (22), we can also obtain the same dispersion relation in eq (25a) for a general anisotropic medium. But there will be some errors when using the rotational waveform directly measured by the sensor. An alternative method is to use the array derived rotation (ADR, Spudich et al. 1995; Spudich & Fletcher 2008) approach to retrieve the rotational waveform. However, it requires that the stations should be distributed as uniformly as possible in different azimuths and its frequency range is also limited by the aperture of a seismic array. It does not need to calculate the vertical partial derivative $\partial_z$, which means the $V'$ in eq. (28) is equal to zero. Then we can also obtain the same dispersion eqs (25a) for accurately calculating the velocity of qR wave by combining eqs (28) and (22). This method has been verified in our recent real data case (Tang et al. 2023).

In this paper, we only study the fundamental mode due to the following two reasons: (1) the perturbation theory used to derive the dispersion formulas is valid only for the fundamental mode at short periods, and for fundamental modes and lowest overtones at longer periods (Maupin 1989). The Rayleigh wave fundamental mode has a dispersion curve isolated enough from those of the other modes not to suffer from strong coupling with neighbouring modes, while the Love wave fundamental mode will be affected by fundamental and higher modes of the Rayleigh wave, and the higher modes of Rayleigh wave and Love wave will suffer from coupling with neighbouring modes seriously and (2) another point to be considered is that the dispersion curve can be accurately obtained by using our approach only when there is the energy of single-mode (fundamental mode) in the seismogram. When the higher mode energy of the surface wave is relatively strong, it is necessary to separate the modes before calculating the dispersion, which is described in detail by Tang & Fang (2023) and Kurrle et al. (2010), but it is currently difficult to use a single station to achieve the mode separation. This is the reason why this method is more suitable for long periods and teleseismic records of shallow seismic sources, when the higher modes are much weaker than the fundamental mode that can be neglected. In the future, with enough rotation sensors, we can use the array method (Tang & Fang 2023) to separate and calculate the local dispersion of higher modes.

## 5 CONCLUSION

We derived the expressions for calculating the dispersion curves in weakly anisotropic media, using a single observation point which consists of translational displacement and rotation or strain. We analysed the influence of surface wave coupling, with effects depending on the anisotropic strength. When the radiation amplitude of the wavefield is very small, the coupled wave will seriously affect the measurements of dispersion curves. Therefore, it is necessary to select a waveform with strong radiation amplitude or use multiple seismograms to fit a stable dispersion value. Even in general anisotropic media, the polarization plane of rotation is perpendicular to the propagation direction, which provides a new method to extract the azimuth of wave propagation. We numerically demonstrate the effectiveness and applicability of the frequency-dependent amplitude ratio method for deriving anisotropy in the Earth.

## ACKNOWLEDGMENTS

This work is funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 955515. We thank Prof Xinding Fang for providing the numerical simulation 3-D FD anisotropic code. We thank Prof Valerie Maupin, Prof Toshiro Tanimoto, Prof Huajian Yao and one anonymous reviewer for their constructive comments that help us improve the paper.
DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

REFERENCES

Aki, K. & Richards, P.G., 2002. *Quantitative Seismology*, 2nd edn, University Science Books.

Al-Lazki, A.I., Sandvol, E., Seber, D., Barazangi, M., Turkelli, N. & Mohamad, R., 2004. Pn tomographic imaging of mantle lid velocity and anisotropy at the junction of the Arabian, Eurasian and African plates, *J. geophys. Int.*, **158**(3), 1024–1040.

Anderson, D.L., 1965. Recent evidence concerning the structure and composition of the Earth's mantle, *Phys. Chem. Earth*, **6**, 1–131.

Anderson, D.L. & Dziewonski, A.M., 1982. Upper mantle anisotropy: evidence from free oscillations, *J. geophys. Int.*, **69**(2), 383–404.

Bernauer, F. et al., 2021. Rotation, strain, and translation sensors performance tests with active seismic sources, *Sensors*, **21**(1), 264, doi:10.3390/s21010264.

Botev, Z.I., Grotowski, J.F. & Kroese, D.P., 2010. Kernel density estimation via diffusion, *Ann. Stat.*, **38**(5), 2916–2957.

Brown, M., Lin, G., Matsuzawa, H. & Yoshizawa, K., 2022. Recovery of love wave overtone waveforms and dispersion curves from single-station seismograms using time-warping, *J. geophys. Int.*, **230**(1), 70–83.

Chen, X., 1993. A systematic and efficient method of computing normal modes for multilayered half-space, *J. geophys. Int.*, **115**(2), 391–409.

Chen, X., 1999. Seismogram synthesis in multi-layered half-space, *Earthq. Res. China*, **13**(2), 150–174.

Chou, P.C. & Pagano, N.J., 1992. *Elasticity: Tensor, Dyadic, and Engineering Approaches*, Courier Corporation.

Crampin, S., 1977. A review of the effects of anisotropic layering on the propagation of seismic waves, *J. geophys. Int.*, **49**(1), 9–27.

Crampin, S. & Chastin, S., 2003. A review of shear wave splitting in the crack-critical crust, *J. geophys. Int.*, **155**(1), 221–240.

Ekström, G., Tromp, J. & Larson, E.W., 1997. Measurements and global models of surface wave propagation, *J. geophys. Res.*, **102**(B4), 8137–8157.

Fang, X. & Tang, L., 2021. A single station approach for subsurface anomaly detection, in *SEG/AAPG/SEPM First International Meeting for Applied Geoscience & Energy*, OnePetro.

Fang, X., Fehler, M.C., Zhu, Z., Zheng, Y. & Burns, D.R., 2014. Reservoir fracture characterization from seismic scattered waves, *J. geophys. Int.*, **196**(1), 481–492.

Fichtner, A. & Igel, H., 2009. Sensitivity densities for rotational ground-motion measurements, *Bull. seism. Soc. Am.*, **99**(2B), 1302–1314.

Forsyth, D.W., 1975. The early structural evolution and anisotropy of the oceanic upper mantle, *J. geophys. Int.*, **43**(1), 103–162.

Forsyth, D.W., Li, A., Levander, A. & Nolet, G., 2005. Array analysis of two-dimensional variations in surface wave phase velocity and azimuthal anisotropy in the presence of multipathing interference, in *Seismic Earth: Array Analysis of Broadband Seismograms, Geophys. Monogr. Ser.*, Vol. 157, pp. 81–97, eds Levander, A. & Nolet, G., AGU.

Francis, T., 1969. Generation of seismic anisotropy in the upper mantle along the mid-oceanic ridges, *Nature*, **221**(5176), 162–165.

Gupta, I.N., 1973. Premonitory variations in S-wave velocity anisotropy before earthquakes in Nevada, *Science*, **182**(4117), 1129–1132.

Hess, H., 1964. Seismic anisotropy of the uppermost mantle under oceans, *Nature*, **203**(4945), 629–631.

Igel, H., Cochard, A., Wassermann, J., Flaws, A., Schreiber, U., Velikosel'tsev, A. & Pham Dinh, N., 2007. Broad-band observations of earthquake-induced rotational ground motions, *J. geophys. Int.*, **168**(1), 182–196.

Igel, H. et al., 2021. ROMY: a multicomponent ring laser for geodesy and geophysics, *J. geophys. Int.*, **225**(1), 684–698.

Keil, S., Wassermann, J. & Igel, H., 2021. Single-station seismic microzonation using 6C measurements, *J. Seismol.* **25**(1), 103–114.

Kurrle, D., Igel, H., Ferreira, A.M., Wassermann, J. & Schreiber, U., 2010. Can we estimate local love wave dispersion properties from collocated

amplitude measurements of translations and rotations?, *Geophys. Res. Lett.*, **37**(4), L04307, doi:10.1029/2009GL042215.

Legendre, C.P., Zhao, L. & Tseng, T.-L., 2021. Large-scale variation in seismic anisotropy in the crust and upper mantle beneath Anatolia, Turkey, *Commun. Earth Environ.*, **2**(1), 1–7.

Levshin, A.L., Pisarenko, V. & Pogrebinsky, G., 1972. On a frequency-time analysis of oscillations, *Annales de geophysique*, **28**, 211–218.

Maupin, V., 1989. Surface waves in weakly anisotropic structures: on the use of ordinary or quasi-degenerate perturbation methods, *J. geophys. Int.*, **98**(3), 553–563.

Maupin, V., 2004. Comment on 'the azimuthal dependence of surface wave polarization in a slightly anisotropic medium' by T. Tanimoto, *J. geophys. Int.*, **159**(1), 365–368.

Mochizuki, E., 1986. The free oscillations of an anisotropic and heterogeneous Earth, *J. geophys. Int.*, **86**(1), 167–176.

Montagner, J.-P. & Nataf, H.-C., 1986. A simple method for inverting the azimuthal anisotropy of surface waves, *J. geophys. Res.*, **91**(B1), 511–520.

Montagner, J.-P. & Tanimoto, T., 1991. Global upper mantle tomography of seismic velocities and anisotropies, *J. geophys. Res.*, **96**(B12), 20 337–20 351.

Mutlu, A.K. & Karabulut, H., 2011. Anisotropic pn tomography of Turkey and adjacent regions, *J. geophys. Int.*, **187**(3), 1743–1758.

Nataf, H.-C., Nakanishi, I. & Anderson, D.L., 1984. Anisotropy and shear-velocity heterogeneities in the upper mantle, *Geophys. Res. Lett.*, **11**(2), 109–112.

Noe, S., Yuan, S., Montagner, J. & Igel, H., 2022. Anisotropic elastic parameter estimation from multicomponent ground-motion observations: a theoretical study, *J. geophys. Int.*, **229**(2), 1462–1473.

Pancha, A., Webb, T., Stedman, G., McLeod, D. & Schreiber, K., 2000. Ring laser detection of rotations from teleseismic waves, *Geophys. Res. Lett.*, **27**(21), 3553–3556.

Park, J. & Yu, Y., 1993. Seismic determination of elastic anisotropy and mantle flow, *Science*, **261**(5125), 1159–1162.

Pedersen, H.A., 2006. Impacts of non-plane waves on two-station measurements of phase velocities, *J. geophys. Int.*, **165**(1), 279–287.

Schmelzbach, C. et al., 2018. Advances in 6C seismology: applications of combined translational and rotational motion measurements in global and exploration seismology, *Geophysics*, **83**(3), WC53–WC69.

Schreiber, U., Igel, H., Cochard, A. & Velikoseltsev, A., 2005. The GEOsensor project: rotations — a new observable for seismology, in *Observation of the Earth System from Space*, eds Flurry, J., Rummel, R., Reigber, C., Rothacher, M., Boedecker, G. & Schreiber, U., Springer, 427pp.

Schulte-Pelkum, V., Masters, G. & Shearer, P.M., 2001. Upper mantle anisotropy from long-period P polarization, *J. geophys. Res.*, **106**(B10), 21 917–21 934.

Shapiro, N.M. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.*, **31**(7), L07614, doi:10.1029/2004GL019491.

Shapiro, N.M., Campillo, M., Stehly, L. & Ritzwoller, M.H., 2005. High-resolution surface-wave tomography from ambient seismic noise, *Science*, **307**(5715), 1615–1618.

Silver, P.G., 1996. Seismic anisotropy beneath the continents: probing the depths of geology, *Ann. Rev. Earth planet. Sci.*, **24**(1), 385–432.

Smith, M.L. & Dahlen, F., 1973. The azimuthal dependence of love and Rayleigh wave propagation in a slightly anisotropic medium, *J. geophys. Res.* **78**(17), 3321–3333.

Sollberger, D., Greenhalgh, S.A., Schmelzbach, C., Van Renterghem, C. & Robertsson, J.O., 2018. 6-C polarization analysis using point measurements of translational and rotational ground-motion: theory and applications, *J. geophys. Int.*, **213**(1), 77–97.

Spudich, P. & Fletcher, J.B., 2008. Observation and prediction of dynamic ground strains, tilts, and torsions caused by the Mw 6.0 2004 Parkfield,
California, earthquake and aftershocks, derived from UPSAR array observations, *Bull. seism. Soc. Am.*, **98**(4), 1898–1914.

Spudich, P., Steck, L.K., Hellweg, M., Fletcher, J. & Baker, L.M., 1995. Transient stresses at Parkfield, California, produced by the *M* 7.4 Landers earthquake of June 28, 1992: observations from the UPSAR dense seismograph array, *J. geophys. Res.*, **100**(B1), 675–690.

Stehly, L., Fry, B., Campillo, M., Shapiro, N., Guilbert, J., Boschi, L. & Giardini, D., 2009. Tomography of the alpine region from observations of seismic ambient noise, *J. geophys. Int.*, **178**(1), 338–350.

Takeuchi, H. & Saito, M., 1972. Seismic surface waves, *Methods Comput. Phys.*, **11**, 217–295.

Tang, L. & Fang, X., 2021a. Generation of 6C synthetic seismograms in stratified vertically transversely isotropic media using a generalized reflection and transmission coefficient method, *J. geophys. Int.*, **225**(3), 1554–1585.

Tang, L. & Fang, X., 2021b. Numerical study of the application of six-component ambient seismic noise data for high-resolution imaging of lateral heterogeneities, in *SEG/AAPG/SEPM First International Meeting for Applied Geoscience & Energy*, OnePetro.

Tang, L. & Fang, X., 2023. Application of six-component ambient seismic noise data for high-resolution imaging of lateral heterogeneities, *J. geophys. Int.*, **232**(3), 1756–1784.

Tang, L., Igel, H. & Montagner, J.-P., 2023. Seismic anisotropy from 6C observations, Tech. rep., Copernicus Meetings.

Tanimoto, T., 1986. Free oscillations of a slightly anisotropic Earth, *J. geophys. Int.*, **87**(2), 493–517.

Tanimoto, T., 2004. The azimuthal dependence of surface wave polarization in a slightly anisotropic medium, *J. geophys. Int.*, **156**(1), 73–78.

Tanimoto, T. & Anderson, D.L., 1985. Lateral heterogeneity and azimuthal anisotropy of the upper mantle: Love and Rayleigh waves 100–250 s, *J. geophys. Res.*, **90**(B2), 1842–1858.

Teanby, N., Kendall, J.-M., Jones, R. & Barkved, O., 2004. Stress-induced temporal variations in seismic anisotropy observed in microseismic data, *J. geophys. Int.*, **156**(3), 459–466.

Vinnik, L.P., Farra, V. & Romanowicz, B., 1989. Azimuthal anisotropy in the earth from observations of SKS at GEOSCOPE and NARS broadband stations, *Bull. seism. Soc. Am.*, **79**(5), 1542–1558.

Wassermann, J., Wietek, A., Hadzioinnou, C. & Igel, H., 2016. Toward a single-station approach for microzonation: using vertical rotation rate to estimate Love-wave dispersion curves and direction finding, *Bull. seism. Soc. Am.*, **106**(3), 1316–1330.

Wassermann, J., Braun, T., Ripepe, M., Bernauer, F., Guattari, F. & Igel, H., 2022. The use of 6DOF measurement in volcano seismology – a first application to Stromboli Volcano, *J. Volc. Geotherm. Res.*, **424**, 107499, doi:10.1016/j.jvolgeores.2022.107499.

Yao, H., van Der Hilst, R.D. & De Hoop, M.V., 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—I. Phase velocity maps, *J. geophys. Int.*, **166**(2), 732–744.

Yao, H., Beghein, C. & Van Der Hilst, R.D., 2008. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis—II. Crustal and upper-mantle structure, *J. geophys. Int.*, **173**(1), 205–219.

Yao, H., Van Der Hilst, R.D. & Montagner, J.-P., 2010. Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography, *J. geophys. Res.*, **115**(B12), B12307, doi:10.1029/2009JB007142.

Yuan, S., Gessele, K., Gabriel, A.-A., May, D.A., Wassermann, J. & Igel, H., 2021. Seismic source tracking with six degree-of-freedom ground motion observations, *J. geophys. Res.*, **126**(3), e2020JB021112.

Zembaty, Z., Bernauer, F., Igel, H. & Schreiber, K.U., 2021. Rotation rate sensors and their applications, *Sensors*, **21**(16), 5344, doi:10.3390/s21165344.

## APPENDIX A

The elastic tensor of VTI media can be expressed as:

$$C_{VTI} = \begin{pmatrix} C_{11} & C_{11} - 2 \cdot C_{66} & C_{13} & 0 & 0 & 0 \\ C_{11} & C_{13} & 0 & 0 & 0 & 0 \\ C_{33} & 0 & 0 & 0 & 0 & 0 \\ & C_{44} & 0 & 0 & 0 & 0 \\ \text{sym.} & & C_{44} & 0 & 0 & 0 \\ & & & C_{66} & 0 & 0 \end{pmatrix} \quad (A1)$$

The elastic tensor of orthorhombic media can be expressed as:

$$C_{\text{orth}} = \begin{pmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{22} & C_{23} & 0 & 0 & 0 & 0 \\ C_{33} & 0 & 0 & 0 & 0 & 0 \\ & C_{44} & 0 & 0 & 0 & 0 \\ \text{sym.} & & C_{55} & 0 & 0 & 0 \\ & & & C_{66} & 0 & 0 \end{pmatrix}. \quad (A2)$$

The elastic tensor of monoclinic media can be expressed as:

$$C_{\text{mono}} = \begin{pmatrix} C_{11} & C_{12} & C_{13} & 0 & 0 & C_{16} \\ C_{22} & C_{23} & 0 & 0 & 0 & C_{26} \\ C_{33} & 0 & 0 & 0 & 0 & C_{36} \\ & C_{44} & C_{45} & 0 & 0 & 0 \\ \text{sym.} & & C_{55} & 0 & 0 & 0 \\ & & & C_{66} & 0 & 0 \end{pmatrix}. \quad (A3)$$

## APPENDIX B

In this appendix, we introduce the linear regression least-square solution of the ratio equation.

We first use the time-frequency transformation method (e.g. the wavelet transformation) to obtain the time-frequency spectra of translational displacement and rotation or radial strain, which have been successfully used in previous studies (Igel et al. 2007; Kurrle et al.
2010), and then combine the least-square algorithm to calculate the dispersion curves. We use the dispersion eq. (17) of isotropic media as an example to demonstrate how to calculate a stable phase velocity. In the time-frequency spectra, eq. (17) can be expressed as

$$c_L(\omega) = \left| \frac{\dot{u}_t(\omega, t_i)}{2\Omega_z(\omega, t_i)} \right|, \qquad (B1)$$

where $t_i(i = 1...N)$ represents any time point in the time window for qL wave. $N$ is the number of data points. Eq. (B1) is written in the form of a matrix

$$G_0 m = d_0, \qquad (B2)$$

where $G_0 = 2[|\Omega_z(\omega, t_1)|, ..., |\Omega_z(\omega, t_N)|]^T$, $m = c_L(\omega)$ and $d_0 = [|\dot{u}_t(\omega, t_1)|, ..., |\dot{u}_t(\omega, t_N)|]^T$. To select the appropriate data point, we introduce a weight function $wf(\omega, t_i)$ for each data point whose value is equal to 1 or 0. Consequently, eq. (B2) can be given as

$$G_m = d, \qquad (B3)$$

where $G = 2[wf(\omega, t_1)|\Omega_z(\omega, t_1)|, ..., wf(\omega, t_N)|\Omega_z(\omega, t_N)|]^T$ and $d = [wf(\omega, t_1)|\dot{u}_t(\omega, t_1)|, ..., wf(\omega, t_N)|\dot{u}_t(\omega, t_N)|]^T$. Consequently, the least-square solution of eq. (B3) can be expressed as;

$$c_L(\omega) = \frac{\sum_{i=1}^{N} wf(\omega, t_i)^2 |\Omega_z(\omega, t_i) \dot{u}_t(\omega, t_i)|}{2 \sum_{i=1}^{N} wf(\omega, t_i)^2 |\Omega_z(\omega, t_i)|^2}, \qquad (B4)$$

when $wf(\omega, t_i)$ is equal to 0 which means that we remove this data point.

## APPENDIX C

In isotropic media, the polarization of rotation and strain at the surface ($z = 0$) can be expressed in Cartesian coordinate as:

For the Love wave;

$$\left\{ \begin{array}{l} \Omega_x = 0 \\ \Omega_y = 0 \\ \Omega_z = i\kappa_{l0}W/2 \end{array} \right. \qquad (C1)$$

$$\left\{ \begin{array}{l} \varepsilon_{xx} = -i\kappa_{l0} \sin \psi \cos \psi W \\ \varepsilon_{yy} = i\kappa_{l0} \sin \psi \cos \psi W \\ \varepsilon_{zz} = 0 \\ \varepsilon_{xy} = iW\kappa_{l0}(\cos \psi^2 - \sin \psi^2)/2 \\ \varepsilon_{xz} = 0 \\ \varepsilon_{yz} = 0, \end{array} \right. \qquad (C2)$$

For the Rayleigh wave;

$$\left\{ \begin{array}{l} \Omega_x = -\kappa_{r0} \sin \psi U \\ \Omega_y = \kappa_{r0} \cos \psi U \\ \Omega_z = 0 \end{array} \right. \qquad (C3)$$

$$\left\{ \begin{array}{l} \varepsilon_{xx} = i\kappa_{r0} \cos \psi^2 V \\ \varepsilon_{yy} = i\kappa_{r0} \sin \psi^2 V \\ \varepsilon_{zz} = -\frac{\lambda}{(\lambda+2\mu)} i\kappa_{r0} V \\ \varepsilon_{xy} = i\kappa_{r0} \sin \psi \cos \psi V \\ \varepsilon_{xz} = 0 \\ \varepsilon_{yz} = 0, \end{array} \right. \qquad (C4)$$

where we again omit the exponential term $e^{ik(\cos\psi x + \sin\psi y) - i\omega t}$. For the sake of simplicity, the following derivation will also omit the term.

In orthorhombic media, the polarization of rotation and strain at the surface ($z = 0$) can be expressed in Cartesian coordinate as

For the qL wave;

$$\left\{ \begin{array}{l} \Omega_x = -\kappa_l \sin \psi T U \\ \Omega_y = \kappa_l \cos \psi T U \\ \Omega_z = i\kappa_l W/2 \end{array} \right. \qquad (C5)$$

$$\left\{ \begin{array}{l} \varepsilon_{xx} = -i\kappa_l \cos \psi (\cos \psi T V - \sin \psi W) \\ \varepsilon_{yy} = i\kappa_l \sin \psi (\sin \psi T V + \cos \psi W) \\ \varepsilon_{zz} = -\frac{C_{33}^2}{C_{33}} i\kappa_l \cos \psi (\cos \psi T V - \sin \psi W) - \frac{C_{33}^2}{C_{33}} i\kappa_l \sin \psi (\sin \psi T V + \cos \psi W) \\ \varepsilon_{xy} = i\kappa_l (2 \sin \psi \cos \psi T V + \cos \psi^2 W - \sin \psi^2 W)/2 \\ \varepsilon_{xz} = 0 \\ \varepsilon_{yz} = 0. \end{array} \right. \qquad (C6)$$
For the qR wave;

$$ 
\left\{
\begin{array}{l}
\Omega_x = -\kappa_r \sin \psi U \\
\Omega_y = \kappa_r \cos \psi U \\
\Omega_z = -i\kappa_r T W/2
\end{array}
\right.
\qquad (C7) $$

$$ 
\left\{
\begin{array}{l}
\varepsilon_{xx} = i\kappa_r \cos \psi(\cos \psi V + \sin \psi T W) \\
\varepsilon_{yy} = i\kappa_r \sin \psi(\sin \psi V - \cos \psi T W) \\
\varepsilon_{zz} = -\frac{C_{13}}{C_{33}}i\kappa_r \cos \psi(\cos \psi V + \sin \psi T W) - \frac{C_{23}}{C_{33}}i\kappa_r \sin \psi(\sin \psi V - \cos \psi T W) \\
\varepsilon_{xy} = i\kappa_r(2 \sin \psi \cos \psi V - \cos \psi^2 T W + \sin \psi^2 T W)/2 \\
\varepsilon_{xz} = 0 \\
\varepsilon_{yz} = 0.
\end{array}
\right.
\qquad (C8) $$

In general anisotropic media, the polarization of rotation and strain at the surface ($z=0$) can be expressed in Cartesian coordinate as
For the qL wave;

$$ 
\left\{
\begin{array}{l}
\Omega_x = -(\kappa_l \sin \psi T U + \cos \psi W' + \sin \psi T V')/2 \\
\Omega_y = (\cos \psi T V' + \kappa_l \cos \psi T U - \sin \psi W')/2 \\
\Omega_z = i\kappa_l W/2
\end{array}
\right.
\qquad (C9) $$

$$ 
\left\{
\begin{array}{l}
\varepsilon_{xx} = i\kappa_l \cos \psi(\cos \psi T V - \sin \psi W) \\
\varepsilon_{yy} = i\kappa_l \sin \psi(\sin \psi T V + \cos \psi W) \\
\varepsilon_{zz} = -\frac{i\kappa_l}{C_{33}}[C_{13}(-\sin \psi \cos \psi W + \cos \psi^2 T V) + C_{23}(\sin \psi \cos \psi W + \sin \psi^2 T V) + \frac{C_{36}}{2}(2 \sin \psi \cos \psi T V + \cos \psi^2 W - \sin \psi^2 W)] \\
\varepsilon_{xy} = i\kappa_l(2 \sin \psi \cos \psi T V + \cos \psi^2 W - \sin \psi^2 W)/2 \\
\varepsilon_{xz} = (\cos \psi T V' - \kappa_l \cos \psi T U - \sin \psi W')/2 \\
\varepsilon_{yz} = (\cos \psi W' + \sin \psi T V' - T\kappa_l \sin \psi U)/2.
\end{array}
\right.
\qquad (C10) $$

For the qR wave;

$$ 
\left\{
\begin{array}{l}
\Omega_x = -(\kappa_r \sin \psi U - \cos \psi T W' + \sin \psi V')/2 \\
\Omega_y = (\cos \psi V' + \kappa_r \cos \psi U + \sin \psi T W')/2 \\
\Omega_z = -i\kappa_r T W/2
\end{array}
\right.
\qquad (C11) $$

$$ 
\left\{
\begin{array}{l}
\varepsilon_{xx} = i\kappa_r \cos \psi(\cos \psi V + \sin \psi T W) \\
\varepsilon_{yy} = i\kappa_r \sin \psi(\sin \psi V - \cos \psi T W) \\
\varepsilon_{zz} = -\frac{i\kappa_r}{C_{33}}[C_{13}(\cos \psi^2 V + \sin \psi \cos \psi T W) + C_{23}(-\sin \psi \cos \psi T W + \sin \psi^2 V) + \frac{C_{36}}{2}(2 \sin \psi \cos \psi V + (\sin \psi^2 - \cos \psi^2) T W)] \\
\varepsilon_{xy} = i\kappa_r(2 \sin \psi \cos \psi V - \cos \psi^2 T W + \sin \psi^2 T W)/2 \\
\varepsilon_{xz} = (\cos \psi V' - \kappa_r \cos \psi U + \sin \psi T W')/2 \\
\varepsilon_{yz} = (-\cos \psi T W' + \sin \psi V' - \kappa_r \sin \psi U)/2.
\end{array}
\right.
\qquad (C12) $$

It should be noted that under the free surface boundary condition, $W'(0) = \partial_z W|_{z=0} = \sigma_{yz} = 0$ in the reference isotropic media (Takeuchi & Saito 1972).