Artifacts removal from multimode dispersion image: new imaging conditions of beamforming and random array

SUMMARY

Beamforming (BF) has been demonstrated to extract multimode surface wave dispersion curves from ambient seismic noise. However, due to the limited sampling of the array and the complex distribution of the noise sources, the dispersion image generated by the array-based technique is usually contaminated by aliasing or artifacts. According to seismic interferometry theory, the Green's function (GF) in the time domain can be retrieved using the noise cross-correlation function (NCF). The Fourier transform of NCFs, that is, the spatial coherence function, is related to the imaginary part of the frequency domain GF. For the vertical component of the surface wave, it corresponds to the zero-order Bessel function of the first kind, that is, the standing wave containing propagating waves in two directions described by positive and negative vector wavenumber. In array techniques based on wavefield transforms, it is common to adopt the propagating wave instead of the standing wave to eliminate the aliasing associated with the negative wavenumber, that is, to replace the Bessel function using the Hankel function or to construct a complete GF via the Hilbert transform. In this paper, we quantitatively analyse the characteristics of three types of aliasing, that is, the aliasing associated with the period extension of the positive wavenumber, the aliasing associated with the negative wavenumber and those associated with the constant wavenumber. The theoretical representations of different imaging conditions are derived for the finite sampling of the wavefield. A new BF imaging condition is then proposed to remove the crossed artifacts, a type of aliasing associated with the negative wavenumber. The new imaging condition relies only on the computed NCFs and does not require reconstruction of the complete GF via the Hilbert transform. The advantage of random sampling in removing artifacts is illustrated. A random array design scheme is suggested by investigating the array performance of the random array and the array designed using tiles of the Hat family newly discovered in the field of monotile aperiodic tiling. We show the artifacts associated with the constant wavenumber, which are usually manifested as a straight line in the dispersion image of the frequency–velocity domain, also known as radial artifacts, can be eliminated by windowing the NCFs.

Image processing, Seismic interferometry, Seismic noise, Surface waves and free oscillations

1 INTRODUCTION

For decades, seismic interferometry (SI; e.g. Nakata et al. 2019) and dense array observation have been two important advances in imaging technology and observation in seismology. This makes it possible to apply an array-based method to subsets of the dense array and extract the phase velocity of the surface waves from ambient seismic noise, and thereby to directly obtain the lateral variation of the phase velocity, without tomography (e.g. Wang et al. 2020; Li et al. 2021; Fu et al. 2022; Qin et al. 2022; Liu et al. 2023; Tsarsitalidou et al. 2024). In particular, it has been illustrated that the higher normal modes or even leaky modes can be extracted by array-based techniques such as wavefield transform (Fourier–Bessel transform or F-J) (e.g. Forbriger 2003; Wang et al. 2019; Hu et al. 2020), Beamforming (BF; e.g. Qin et al. 2022; Wu et al. 2023b), and Spatial Autocorrelation (e.g. Aki 1957; Yamaya et al. 2021). The array-based imaging method of surface waves is becoming increasingly popular (e.g. Yamaya et al. 2021; Chen et al. 2023; Wu et al. 2023a, 2023b) and is changing the dominance of the two-station surface wave method. Since the multimode dispersion curves or the energy spectrum is crucial for the surface wave inversion, how to generate the multimode dispersion image with high resolution from ambient seismic noise has attracted the attention of researchers in seismological community (e.g. Cheng et al. 2021; Takagi & Nishida 2022; Nimiya et al. 2023; Viens et al. 2023; Yao et al. 2023; Qin & Lu 2024).

The stations involved in the array is usually finite in the application. In fact, it is almost impossible and in general unnecessary to sample the data fully in time and space to reduce the burden of data acquisition and storage, especially for the massive continuous seismic noise recordings. The finite sampling on the wavefield contradicts the assumption of full collection required by most signal-processing techniques (Snieder & Wakin 2022). The problem posed by insufficient sampling is that the dispersion image generated by the array-based technique has some artifacts due to aliasing, which affects the identification of the actual mode branches of the dispersion curves.

One approach to eliminate these artifacts is to predict them theoretically by investigating their causes, and then mitigate or eliminate them from the dispersion image by modifying the imaging conditions. For example, when processing the seismogram gathered from active source, Forbriger (2003) replaced the Bessel function with Hankel function in the F-J transform to eliminate aliasing arising from negative wavenumber representing inward-propagating waves. Xi et al. (2021) and Zhou & Chen (2022) generalized this principle to virtual records represented by noise cross-correlation functions (NCFs). Since NCF corresponds to the imaginary part of Green's function (GF) in the frequency domain (e.g. Lobkis & Weaver 2001; Campillo & Paul 2003; Wapenaar & Fokkema 2006; Lu 2021), they construct the complete GF (real and imaginary parts) by Hilbert transform of NCFs (e.g. Nakahara 2006), and then modify the F-J transform using the complete GF to remove the crossed artifacts.

In this paper, a new BF imaging condition is proposed to eliminate artifacts caused by aliasing. The new condition utilizes the real and imaginary parts of NCFs. The construction of the complete GF is not required. BF is an array-based technique for obtaining multimode surface waves from ambient seismic noise. Different from the numerical integration in F-J, delay and sum are performed in BF to match the slowness and direction of an incident plane wave (e.g. Gal, M. & Reading 2019). The azimuth-averaged or azimuth-dependent phase velocities of the wave propagating in the structure can be estimated using the slowness and azimuth associated with the maximum of the beamformer output (e.g. Harmon et al. 2008; Roux & Ben-Zion 2017; Löer et al. 2018; Wang et al. 2020; Qin et al. 2022). Some improved schemes are proposed to enhance the resolution of the dispersion image and reduce the effect of the irregular array geometry and anisotropic sources (Qin & Lu 2024). However, artifacts also appear in the dispersion images generated by BF, especially for the regular array geometry with equal spacing. To remove artifacts, we first investigate quantitatively three types of aliasing, that is, the aliasing associated with the periodic extension of the positive wavenumber, the aliasing associated with the periodic extension of the negative wavenumber, and those associated with the constant wavenumber. The theoretical representations of different imaging conditions are then derived for the finite sampling of the wavefield. A new imaging condition is then proposed to eliminate artifacts associated with the negative wavenumber. The new imaging condition relies only on the computed NCFs and does not require reconstruction of the complete GF via the Hilbert transform.

According to the SI theory, GF between stations A and B can be approximated by NCF. The NCF contains causal and acausal parts, which corresponds to waves travelling from A to B, or from B to A, respectively. For the 2-D case, which is common for noise-based surface wave method, the real part of the causal NCF corresponds to |${J_0}(kr)$|⁠, the zero-order Bessel function of the first kind, while its imaginary part corresponds to the zero-order Struve function |${S_0}(kr)$| (e.g. Tsai 2009; Boschi et al. 2013; Lu 2021). The proposed new imaging conditions are constructed directly based on the real and imaginary parts of the causal NCFs.

Another way to reduce artifacts is to replace regular sampling with random sampling since aliasing associated with periodic sampling would be avoided for random sampling (Snieder & Wakin 2022). Inspired by advances in the field of aperiodic tiling (Smith et al. 2024), Mordret & Grushin (2024) investigated the performance of seismic array designed by the tiles of the newly discovered Hat family, assuming that the stations are located on the vertices of tiles. They show that the monotile aperiodic seismic arrays outperform regular arrays, while Specter tiles performs better than the other tiles of the Hat family, such as the Hat and the Turtle, etc. Monotile aperiodic seismic array may provide a solution in seismic array design. However, as they mentioned, ‘the aperiodic nature of the Specter tiling is not the fundamental reason it performs so well’. Aperiodic tiling does not imply non-periodicity in the sampling on spatial wavefield. The array response function (ARF) of the Hat family and several random arrays are investigated in this paper. Some examples of synthetic and field data are presented to show that the array with random station distribution facilitates the elimination of artifacts associated with the periodic extension of the positive and negative wavenumber. As the artifacts associated with the constant wavenumber, which are usually manifested as a straight line, also known as radial artifacts, we show they can be eliminated by pre-processing NCFs.

The paper is structured as follows: the BF schemes for extracting multimode dispersion curves is briefly reviewed in Section 2. In Section 3, the uncertainty of the dispersion measurement and aliasing are investigated based on the ARF. Three types of aliases associated with limited sampling are formularized. The schemes for eliminating crossed artifacts are outlined based on SI theory in Section 4. The new BF imaging conditions are then proposed in this section. In Section 5, the examples of synthetic and field data are provided to illustrate the removal of crossed artifacts using the new imaging condition. In Section 6, the advantage of random sampling in removing artifacts is illustrated. A random array design scheme is then suggested by investigating the array performance of the random array and the array designed by the tiles of Hat family. The removal of radial artifacts is discussed in Section 7. The conclusions and perspectives are given in Section 8.

2 BEAMFORMING

2.1 Beamforming schemes for extracting multimode dispersion curves

Beamforming performs phase correction of the recordings at N stations using a plane wave with a supposed velocity. The beamformer output is calculated by conducting a delay-and-sum of the recordings over the stations. When the supposed phase velocity matches the true one under the array, the beamformer output reaches its maximum and the phase velocity under the array is thereby obtained. The cross-correlation beamforming (CBF) of seismic noise can be written as (e.g. Wang et al. 2020):

$$\begin{eqnarray} \mathrm{ CBF}\left( {k,\omega ,\theta } \right)& =& \frac{1}{{{N^2}}}{\left| {\sum\limits_i^N {u\left( {{{\bf{x}}_i},\omega } \right){e^{ - \iota {\bf{k}} \cdot {{\bf{x}}_i}}}} } \right|^2} \\ &=& \frac{1}{{{N^2}}}{\left( {\sum\limits_i^N {u\left( {{{\bf{x}}_i},\omega } \right){e^{ - \iota {\bf{k}} \cdot {{\bf{x}}_i}}}} } \right)^*}\left( {\sum\limits_j^N {u\left( {{{\bf{x}}_j},\omega } \right){e^{ - \iota {\bf{k}} \cdot {{\bf{x}}_j}}}} } \right) \\ &=& \frac{1}{{{N^2}}}\left| {\sum\limits_i^N {\sum\limits_j^N {{u^*}\left( {{{\bf{x}}_i},\omega } \right)u\left( {{{\bf{x}}_j},\omega } \right){e^{ - \iota {\bf{k}} \cdot ({{\bf{x}}_j} - {{\bf{x}}_i})}}} } } \right| \\ &=& \frac{1}{{{N^2}}}\left| {\sum\limits_i^N {\sum\limits_j^N {{C_{ij}}\left( \omega \right){e^{ - \iota k{r_{ij}}\cos \left( {{\theta _{ij}} - \theta } \right)}}} } } \right| ,\end{eqnarray}$$

(1)

where N is the number of stations in the array. |$u({{\bf{x}}_i},\omega )$| is the Fourier transform of the velocity recordings pertaining to the vertical component at station i, |$k = | {\bf{k}} |$| and |${\bf{k}}$| is the vector wavenumber of the incident plane wave with azimuth |$\theta $|⁠. |$\iota = \sqrt { - 1} $| is the imaginary unit. * represents the complex conjugation. |${r_{ij}} = | {{{\bf{x}}_j} - {{\bf{x}}_i}} |$| represents the interstation distance of the station pair. |${\theta _{ij}}$| represents the azimuth from station i to station j, measured clockwise from the north. The |${C_{ij}}( \omega )$| is called the cross-spectral density matrix (CSDM) (e.g. Gerstoft & Tanimoto 2007; Riahi and Saenger 2014) or cross-covariance matrix (e.g. Capon 1969; Seydoux et al. 2017), which corresponds to NCFs in the time domain and is the element used in traditional noise interferometry. In this paper, it is assumed the records u is the vertical component. Only the cross-correlation of the vertical–vertical (ZZ) components is considered. The ZZ component corresponds to the Rayleigh waves, which propagates as the cylindrical wave and can be approximated by the plane wave in far field. A slight modification can be performed to extend BF for estimation of Love waves using the transverse–transverse (TT).

The factor |$\sqrt {kr} $| is multiplied in weighted cross-correlation beamforming (WCBF) to correct the geometrical spread of NCFs and to improve the resolution of the dispersion image at high frequencies. It is written as (Qin & Lu 2024):

$$\begin{eqnarray} \begin{array}{@{}c@{}} \textrm{WCBF}\left( {k,\omega ,\theta } \right) = \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {{C_{ij}}\left( \omega \right){e^{ - \iota {\bf{k}} \cdot {{\bf{r}}_{ij}}}}} } \left| {\bf{k}} \right|{r_{ij}}\\ = \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {{C_{ij}}\left( \omega \right){e^{ - \iota k{r_{ij}}\cos \left( {{\theta _{ij}} - \theta } \right)}}} } k{r_{ij}} \end{array} . \end{eqnarray}$$

(2)

For CBF and WCBF, the plane wave is projected into all the available interstation. It implies the structure under the array is assumed to be isotropic. If CBF and WCBF are used to estimate the azimuth-dependent phase velocity, the result would be biased due to the array geometry and uneven source distribution (e.g. Lu et al. 2018). A modified beamforming (MCBF) is then proposed to estimate the azimuth-dependent velocity with reducing effect of the irregular geometry and anisotropic sources. The MCBF can be expressed as (Qin & Lu 2024):

$$\begin{eqnarray} \textrm{MCBF}\left( {k,\omega ,\theta = {\theta _{ij}}} \right) = \frac{1}{{N{{({\theta _{ij}})}^2}}}\sum\limits_i^{N({\theta _{ij}})} {\sum\limits_j^{N({\theta _{ij}})} {\sqrt {\pi k{r_{ij}}} {C_{ij}}\left( \omega \right){e^{ - \iota \left( {k{r_{ij}} - \frac{\pi }{4}} \right)}}} } .\end{eqnarray}$$

(3)

In MCBF, for a given interstation orientation |${\theta _{ij}}$|⁠, only the incident plane waves located in the stationary phase region are considered. |$N( {{\theta _{ij}}} )$| represents the number of station pairs located in the stationary phase region. Naturally, MCBF can also be used to estimate the azimuth-averaged velocity. The azimuth-averaged version of MCBF is

$$\begin{eqnarray} \overline {\textrm{MCBF}} (k,\omega ) = \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {\sqrt {\pi k{r_{ij}}} {C_{ij}}(\omega ){e^{ - \iota \left( {k{r_{ij}} - \frac{\pi }{4}} \right)}}} } . \end{eqnarray}$$

(4)

As far as the estimation of the azimuth-averaged velocity, the dispersion image generated by MCBF is similar as that by WCBF.

2.2 Dispersion image generated by different beamforming schemes

As an example, Fig. 1 gives an illustration of the azimuthally averaged dispersion image generated by WCBF and MCBF. The triangles in Fig. 1(a) are the subset of the dense array located in Tongzhou (Beijing) with intervals of 1–2 km (Qin et al. 2022). We first sample the stations so that they are distributed as regularly as possible to highlight the features of the aliasing (Since the random sampling would reduce the structure aliasing, as discussed in Section 6). The 56 stations used in the computation are indicated by blue triangles. The dotted lines in Fig. 1(a) shows the grids with 1 km × 1 km. The used stations are almost located near the nodes. The ZZ component of the NCFs of the interstation are computed using one-month recordings. Figs 1(b) and (c) show the dispersion images generated by WCBF and MCBF, respectively.

Figure 1.

The multimode dispersion images generated by different beamforming schemes. (a) The configuration of the array. The dotted lines in panel (a) shows the grids with 1 km × 1 km. The stations used are the blue triangles consisting of 56 stations, which are sampled as regularly as possible to highlight the artifacts (See Fig. 18 for the results of the grey triangles). (b) The azimuth-averaged dispersion image obtained by WCBF. (c) The dispersion image obtained by MCBF. The typical artifacts predicted by the observed modes based on the eqs (19)–(21) are presented in panel (c). The part enclosed by the solid black line is the artifact related to the constant wavenumber. The orange dash–dotted lines denote crossed artifacts associated with the negative wavenumber. The white dashed (m = 1) and dotted lines (m = −1) denote the artifacts associated with the positive wavenumber for different m. See also Fig. 14(c).

Open in new tab Download slide

The fundamental mode (mode 0) and the first overtone (mode 1) can be observed in Figs 1(b) and (c). The typical artifacts are labelled in Fig. 1(c). These artifacts are predicted based on the observed mode branches (see Section 3.3 for details). Similar to the classic 1-D Fourier transform, the generation of the artifacts in BF is directly related to the discrete sampling of the spatial wavefields, which mainly depends on the array configuration and thereby on the ARF. We first investigate the resolution and artifacts via ARF. A new BF imaging condition is then proposed to reduce the artifacts.

3 ARF, RESOLUTION AND ALIASING

3.1 ARF and resolution

The resolution of BF depends on the array configuration and the characteristics of the wavefield across the array. The wavefield, for example, the energy of the surface wave modes carried by the ambient noise, depends on the structure under the array and source characteristics which are often what we are trying to figure out. Therefore, the resolution we are concerned with and can be improved usually refers to the resolution determined by the array configuration.

The CSDM appearing in eq. (1) is related to the wavefield. Removing CSDM from eq. (1), the remaining part depends on the array configuration only, and can be expressed as (Wathelet et al. 2008; Ruigrok et al. 2017):

$$\begin{eqnarray} \mathrm{ ARF}\left( {\bf{k}} \right) = \frac{1}{{{N^2}}}\left| {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{e^{ - \iota {\bf{k}} \cdot {{\bf{r}}_{ij}}}}} } } \right| . \end{eqnarray}$$

(5)

Since

$$\begin{eqnarray} {e^{ - \iota {\bf{k}} \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}}{\rm{ = }}{\left[ {{e^{ - \iota {\bf{k}} \cdot \left( {{{\bf{x}}_i} - {{\bf{x}}_j}} \right)}}} \right]^ * } . \end{eqnarray}$$

(6)

Eq. (5) can also be written as (Horike 1985; Picozzi et al. 2010):

$$\begin{eqnarray} \mathrm{ ARF}\left( {\bf{k}} \right) = \frac{1}{{{N^2}}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{e^{ - \iota {\bf{k}} \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}}} } . \end{eqnarray}$$

(7)

|$\mathrm{ ARF}({\bf{k}})$| is known as the spatial window function (e.g. Lacoss et al. 1969; Horike 1985) or ARF (e.g. Capon 1969; Rost & Thomas 2002; Ruigrok et al. 2017; Näsholm et al. 2022). Physically it acts as the beamformer output for the wavefield with unit constant amplitude, that is, |${C_{ij}}(\omega ) = 1$| in eq. (1).

For an ideal monochromatic plane wave with angular frequency |${\omega _0}$| and wavenumber |${{\bf{k}}_0}$|⁠, by ignoring the attenuation, the CSDM between stations |${{\bf{x}}_i}$| and |${{\bf{x}}_j}$| can be expressed as (e.g. Löer et al. 2018):

$$\begin{eqnarray} {C_{ij}}\left( \omega \right){\rm{ = }}\delta \left( {\omega {\rm{ - }}{\omega _0}} \right){e^{\iota {{\bf{k}}_0} \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}} , \end{eqnarray}$$

(8)

where |$\delta (\omega - {\omega _0})$| is the Dirac delta function. Substituting eq. (8) into eq. (1) and applying eq. (5), the beamformer output can be written as:

$$\begin{eqnarray} B\left( {\omega ,{\bf{k}}} \right) &=& \frac{1}{{{N^2}}}\left| {\sum\limits_i^N {\sum\limits_j^N {\delta \left( {\omega - {\omega _0}} \right){e^{\iota {{\bf{k}}_0} \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}} \cdot {e^{ - \iota {\bf{k}} \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}}} } } \right| \\ & =& \delta \left( {\omega - {\omega _0}} \right)\frac{1}{{{N^2}}}\left| {\sum\limits_i^N {\sum\limits_j^N {{e^{ - \iota \left( {{\bf{k}} - {{\bf{k}}_0}} \right) \cdot \left( {{{\bf{x}}_j} - {{\bf{x}}_i}} \right)}}} } } \right| \\ & =& \delta \left( {\omega - {\omega _0}} \right)\mathrm{ ARF}\left( {{\bf{k}} - {{\bf{k}}_0}} \right) . \end{eqnarray}$$

(9)

For a more general case with finite length, the estimated wavenumber spectrum is the 2-D convolution of the true spectrum with ARF (corresponding to the product in the frequency domain) (e.g. Lacoss et al. 1969; Asten & Henstridge 1984).

$$\begin{eqnarray} B\left( {\omega ,{\bf{k}}} \right) = {\left| {{W_f}\left( {\omega - {\omega _0}} \right)} \right|^2}\mathrm{ ARF}\left( {{\bf{k}} - {{\bf{k}}_0}} \right) \end{eqnarray}$$

(10)

where |${W_f}(\omega - {\omega _0})$| is the Fourier spectra of the time-series with finite length. The BF resolution is thereby controlled by ARF which depends on the array configuration.

At present, there is no global agreement about the resolving capabilities of an array (Wathelet et al. 2008). For the array with simple and regular geometries, for instance, a linear array with equal spacing, the aliasing and resolution limit can be estimated using the maximum and minimum interstation spacing based on the Nyquist sampling theorem. For the array with irregular geometry, some empirical rules are proposed to estimate the reasonable results achieved by the array. Tokimatsu (1997) use the minimum |${r_{\min }}$| and maximum |${r_{\max }}$| of the interstation distances inside the array to determine the range of the resolved wavelength. As a rule of thumb, the resolved minimum wavelength |${\lambda _{\min }}$| and maximum one |${\lambda _{\max }}$| are respectively |${\lambda _{\min }} = 2{r_{\min }}$| and |${\lambda _{\max }} = 3{r_{\max }}$| (e.g. Tokimatsu 1997; Löer et al. 2020).

A more rigorous definition for resolved wavenumber is based on ARF. For two plane waves travelling across the array with wavenumbers |${{\bf{k}}_1}$| and |${{\bf{k}}_2}$|⁠, the CSDM can be expressed as:

$$\begin{eqnarray} {C_{ij}}(\omega ) = {A_1}(\omega ){e^{\iota {{\bf{k}}_1} \cdot {{\bf{r}}_{ij}}}} + {A_2}(\omega ){e^{\iota {{\bf{k}}_2} \cdot {{\bf{r}}_{ij}}}} . \end{eqnarray}$$

(11)

Where |${A_1}$| and |${A_2}$| represent the amplitudes of two plane waves, respectively. The spectrum in the frequency–wavenumber (f–k) domain estimated from eq. (1) can be written as:

$$\begin{eqnarray} \begin{array}{@{}c@{}} \mathrm{ CBF}({\bf{k}},\omega ,\theta ) = \frac{1}{{{N^2}}}\left| {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {\left( {{A_1}(\omega ){e^{\iota {{\bf{k}}_1} \cdot {{\bf{r}}_{ij}}}} + {A_2}(\omega ){e^{\iota {{\bf{k}}_2} \cdot {{\bf{r}}_{ij}}}}} \right){e^{ - \iota {\bf{k}} \cdot {{\bf{r}}_{ij}}}}} } } \right|\\ \le \frac{1}{{{N^2}}}\left( {\left| {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{A_1}(\omega ){e^{\iota {{\bf{k}}_1} \cdot {{\bf{r}}_{ij}}}}{e^{ - \iota {\bf{k}} \cdot {{\bf{r}}_{ij}}}}} } } \right| + \left| {\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{A_2}(\omega ){e^{\iota {{\bf{k}}_2} \cdot {{\bf{r}}_{ij}}}}{e^{ - \iota {\bf{k}} \cdot {{\bf{r}}_{ij}}}}} } } \right|} \right)\\ = {\left| {{W_1}\left( {\omega - {\omega _1}} \right)} \right|^2}\mathrm{ ARF}\left( {{\bf{k}} - {{\bf{k}}_1}} \right) + {\left| {{W_2}\left( {\omega - {\omega _2}} \right)} \right|^2}\mathrm{ ARF}\left( {{\bf{k}} - {{\bf{k}}_2}} \right) \end{array} .\end{eqnarray}$$

(12)

This means the beamformer output for two plane waves is always lower than the summation of individual plane waves (e.g. Wathelet et al. 2008). The aliasing and resolution of the array can be defined based on the ARF by considering the summation of two shifted ARFs. If two wavenumbers are close to each other, the summation of two shifted ARF would generate a wider main lobe rather than two narrow main lobes. The two wavenumbers would not be distinguished if their difference is less than the width of the main lobe.

Fig. 2(a) shows the distribution of the interstation numbers as a function of spacing and orientation for the array (blue triangles) shown in Fig. 1(a). It shows the samplings of the wavefield are mainly concentrated on some specific points. Fig. 2(b) shows the ARF of the array. The maximum occurred at |${\bf{k}} = 0$|⁠. A series of subpeaks are distributed in the surrounding area. As the wavenumber increases, the amplitude of the subpeaks far away from the main one is gradually weakened. However, the eight secondary peaks have the comparable amplitude with the main one. Fig. 2(c) shows the cross-sections of ARF every two degrees.

Figure 2.

(a) The distribution of the station-pair numbers as a function of spacing and orientation. (b) The ARF of the array in Fig. 1(a) (blue triangles). (c) The cross-sections (grey lines) of ARF every two degrees. The black solid line denotes the azimuth-averaged result. The blue solid line indicates the width of the main lobe at an amplitude of 0.5. The vertical black dashed lines indicate the position of the side lobe.

Open in new tab Download slide

The narrower the central peak of ARF, the more capable the array to distinguish two waves travelling at close wavenumbers (e.g. Wathelet et al. 2008). Different definitions are proposed to evaluate the resolution of BF and the confidence interval that aliasing does not appear. The width of the main lobe, at the edge of which the beamformer output is reduced to a given threshold relative to the maximum, is usually used to define the resolving power of the two wavenumbers. As shown in Fig. 2(c), two wavenumbers cannot be resolved if their difference is less than|$\Delta k$| since the two main lobes would overlap. If the maximum amplitude of the main lobe is 1, the threshold is generally selected as 0.5, that is, the beamformer output reduces to half of the maximum or −3 dB in the logarithmic coordinates (Woods & Lintz 1973; Asten & Henstridge 1984; Wathelet et al. 2008). For this definition, |$\Delta k$| is equal to 0.56 rad km⁻¹ (⁠|$\Delta k$|represented the width indicated by the thick blue line in Fig. 1c) for the array shown in Fig. 1(a).

For an array with a given configuration, the width of the main lobe and thereby the width of the energy belt along the mode branch remains the same for all wavenumbers, and as a result it varies with frequency for a dispersive wave. For the dispersion image in the frequency–velocity (f–v) domain, the upper and lower bounds of the energy belt can be expressed as:

$$\begin{eqnarray} \begin{array}{@{}l@{}} \Delta {v^ + }\left( f \right) = \frac{{{v^2}}}{{\frac{{4\pi f}}{{\Delta k}} - v}}\\ \Delta {v^ - }\left( f \right) = \frac{{{v^2}}}{{\frac{{4\pi f}}{{\Delta k}} + v}} \end{array} . \end{eqnarray}$$

(13)

Where |$\Delta k$| represents the width of ARF main lobe. For the two modes shown in Fig. 1(c), the estimated |$\Delta {v^ + }$| and |$\Delta {v^ - }$| are presented in Fig. 3. As expected, the width of the energy belt is broad at low frequencies and narrow at high. According to eq. (13), both |$\Delta {v^ + }$| and |$\Delta {v^ - }$| tend to zero with increasing frequency. For frequencies greater than a certain value, for example 2.5 Hz for mode 1, the resolving power would disappear due to the absence of the energy belt. Therefore, when generating the dispersion image by BF, the step of the phase velocity in computation should be small enough to ensure that the loss of the resolving power for high frequencies is not caused by the numerical calculation.

Figure 3.

The width of the energy belt varies with frequencies along the mode branches, which is an indication of the accuracy and uncertainty of measurements and can be predicted by eq. (13).

Open in new tab Download slide

According to the sampling theorem, the maximum (⁠|${k_{\max }}$|⁠) and minimum (⁠|${k_{\min }}$|⁠) wavenumber values that can be resolved by the array are can be expressed as:

$$\begin{eqnarray} \begin{array}{@{}l@{}} {k_{\max }} = {\pi {\left/ {\vphantom {\pi {{r_{\min }}}}} \right. } {{r_{\min }}}}\\ {k_{\min }} = {\pi {\left/ {\vphantom {\pi {{r_{\max }}}}} \right. } {{r_{\max }}}} \end{array} , \end{eqnarray}$$

(14)

where |${r_{\max }}$| and |${r_{\min }}$| are the largest and smallest interstation distances, respectively. |${k_{{\rm{max}}}}$| is also called Nyquist wavenumber. For the array shown in Fig. 1a, |${k_{\max }} = 14.1{\rm{ rad\,km^{-1}, }}{k_{{\rm{min}}}} = 0.5{\rm{ rad\,km^{-1}}}$|⁠. However, the range of resolvable wavenumber is not strictly limited to the range given by sampling theorem, particularly for the maximum resolvable wavenumber (e.g. Foti et al. 2018). The range of the resolvable wavenumber also depends on the wavefield and site conditions. In practice, the reliable wavenumber range identified from the dispersion image is usually larger than that given by sampling theorem. For Fig. 3, we empirically estimate that the range of reliable wave numbers is [2k_min, 2k_max]. The resolvable maximum wavenumber is larger than |${k_{\max }}$|⁠, as marked in Fig. 3.

3.2 Subpeaks and aliasing

The |$\mathrm{ ARF}{\rm{(}}{\bf{k}}{\rm{)}}$| is equivalent to an array response to an incident plane wave vertically with a vanished horizontal wavenumber. For a plane wave incident horizontally with wavenumber vector |${{\bf{k}}_0}$|⁠, it can be expressed as |$\mathrm{ ARF}{\rm{(}}{\bf{k}} - {{\bf{k}}_0}{\rm{)}}$|⁠, that is, the spatial shift of |$\mathrm{ ARF}{\rm{(}}{\bf{k}}{\rm{)}}$|⁠. The maximum of the main peak of ARF appears in each azimuth if the sources are uniformly distributed and the orientation of the station-pair has a homogeneous azimuthal coverage. As a result, the beamformer output exhibits an approximately circular energy belt. Fig. 4(a) shows the beamformer output for the array shown in Fig. 1(a) in the wavenumber domain. The circular energy belt centred at the origin corresponds to the main lobe of the ARF with magnitude indicating the relative intensity of the sources (e.g. Gerstoft & Tanimoto 2007). Several circular belts are distributed around the central one. They are associated with the subpeaks of the ARF in Fig. 2(b). At 0.4 Hz, the central energy belt is separated from those around, while they intersect at 0.5 Hz due to insufficient spatial sampling. At 0.7 Hz, two modes exist and hence two circular energy belts are expected in the wavenumber domain. The beamformer output exhibit more complexity due to the cross talk of two belts.

Figure 4.

The beamformer output for the array shown in Fig. 1(a) (blue triangles) at 0.4, 0.5, and 0.7 Hz.

Open in new tab Download slide

Fig. 4(b) shows the results in the phase velocity domain. Since the energy belts associated with aliasing have larger wavenumbers, corresponding to the lower phase velocities, the aliasing associated with them appears in the interior of the main circular belt in the velocity domain, as shown for 0.4 Hz. For 0.5 and 0.7 Hz, the energy belts associated with the subpeaks intersect with the main one of interest, exhibiting more complex characteristics. To generate the dispersion image by averaging the beamformer output over the azimuth, the extremum associated with the main lobe gives the real mode branches, while the extremum associated with the subpeaks is manifested as the aliasing or artifacts.

3.3 Artifacts in the dispersion image

The periodic aliasing in the wavenumber or velocity domain results from the regular sampling (e.g. Zwartjes & Sacchi 2007; Snider & Wakin 2022). At the surface of the isotropic layered model, to sample the cylindrical wavefield using a 2-D array is equivalent to sample them using several regular 1-D arrays with equal spacing but different values. A simple 1-D linear array is hence used to illustrate the artifacts in the dispersion image since it provides a clear picture of artifacts in the dispersion image. The two-layered elastic model shown in Table 1 is considered.

Table 1.

Open in new tab

The parameters of a two-layered model (Zhou & Chen 2022).

Layer thickness (km)	P-wave velocity (km s⁻¹)	S-wave velocity (km s⁻¹)	Density (g cm⁻³)
0.025	1.35	0.2	1.9
\|$\infty $\|	2.0	1	2.5

Table 1.

Open in new tab

The parameters of a two-layered model (Zhou & Chen 2022).

Layer thickness (km)	P-wave velocity (km s⁻¹)	S-wave velocity (km s⁻¹)	Density (g cm⁻³)
0.025	1.35	0.2	1.9
\|$\infty $\|	2.0	1	2.5

As shown in Fig. 5(a), the linear array consists of 12 stations with interval |$\Delta x = 15\,\mathrm{ m}$|⁠. Fig. 5(b) shows the ARF of the array. According to the sampling theorem, the aliasing wavenumber |${k_{\textrm{aliasing}}} = 2\pi /\Delta x = 419\,\mathrm{ rad}\,\mathrm{ km}^{-1}$| appears in the secondary peak. The wavenumber larger than |${k_{\textrm{aliasing}}}$| cannot be resolved and the aliasing would appear at |$m{k_{\textrm{aliasing}}}$|(⁠|$m = 0, \pm 1, \pm 2, \cdots $|⁠) periodically.

$The illustration of the aliasing appears in the dispersion image generated by BF. (a) The array configuration used in the numerical simulation. (b) The ARF. Panels (c) and (d) are the dispersion image in the f–k and f–v domain, respectively. The black dotted lines denote the theoretical dispersion curves. Panels (e) and (f) are the results for a constant CSDM (equivalent to the ARF). The black solid and dashed lines represent $m \cdot {k_{\textrm{aliasing}}}$ and $m/2 \cdot {k_{\textrm{aliasing}}}$, respectively.$

Figure 5.

The illustration of the aliasing appears in the dispersion image generated by BF. (a) The array configuration used in the numerical simulation. (b) The ARF. Panels (c) and (d) are the dispersion image in the f–k and f–v domain, respectively. The black dotted lines denote the theoretical dispersion curves. Panels (e) and (f) are the results for a constant CSDM (equivalent to the ARF). The black solid and dashed lines represent |$m \cdot {k_{\textrm{aliasing}}}$| and |$m/2 \cdot {k_{\textrm{aliasing}}}$|⁠, respectively.

Open in new tab Download slide

The spatial autocorrelation coefficient of the surface wave vertical component can be expressed as (e.g. Haney et al. 2012; Lu 2021):

$$\begin{eqnarray} {C_{ij}}(\omega ) = \frac{1}{{8{c_n}{U_n}I}}{J_0}\left(\frac{\omega }{{{c_n}}}r\right) . \end{eqnarray}$$

(15)

Where |${c_n}$| and |${U_n}$| represent the phase and group velocity of Rayleigh wave with mode n, respectively. I denote the energy integral. In the simulation, the dispersion image can be generated using MCBF by substituting eq. (15) into eq. (4).

The image in the f–k domain presents two types of approximate straight lines that appears periodically, upward-sloping and downward-sloping lines. Upward lines are associated to peaks with positive extremum, while downward lines are characterized by a trough (zero point) between two opposite peaks.

The black dotted lines in Fig. 5(c) denote the mode branches of interest. These rising lines appear periodically as aliasing in the large wave number region in the f–k domain. They behave as artifacts in the f–v domain in the bottom region near zero velocity, as shown in Fig. 5(d). The downward lines in the f–k domain behave as the ‘crossed artifacts’ in the f–v domain. These crossed artifacts intersect with the true mode branch at |$m/2 \cdot {k_{\textrm{aliasing}}}$|⁠, |$m = \pm 1, \pm 2,...$|⁠, as shown by the black dashed lines in Figs 5(c) and (d). Furthermore, in the right-bottom of Fig. 5(c), there are multiple upward lines associated with aliasing, for example, the one appearing from 13 Hz onwards as marked by the black arrow. It is characterized in the f–v domain similarly to the actual mode.

In addition, the artifacts manifested as straight lines in the f–v domain are sometimes observed. These artifacts correspond to a constant wavenumber. As discussed in Section 3.1, if the CSDM is constant, the beamformer output is nothing but the ARF. In the f–k domain, the ARF behaves the horizontal lines positioned at |$m \cdot {k_{\textrm{aliasing}}}$|⁠, |$m = 0, \pm 1, \pm 2,...$|⁠, as shown in Fig. 5(e). These horizontal aliasing will behave as straight lines in the f–v domain, as shown in Fig. 5(f). These aliases are called radial pattern artifacts in Luo et al. (2022). The radial artifacts are not mandatory and can be observed only if the vertically incident wave has the comparable energy magnitude as that of a horizontally incident surface wave.

In summary, for a given angular frequency |$\omega $|⁠, artifacts can be classified into three groups as a function of the actual phase velocity at that frequency. It is assumed the multimode surface wave phase velocity is defined as |${c_n}( \omega )$|⁠, where n represents the mode index. The corresponding wavenumber is |${k_n}( \omega ) = \omega /{c_n}( \omega )$|⁠. The artifacts generated by periodic extension of the actual wavenumber |${k_n}( \omega )$| and the negative wavenumber |$- {k_n}( \omega )$|⁠, as well as those related to the constant wavenumber can be written as:

$$\begin{eqnarray} k_n^{\prime}\left( \omega \right) = {k_n}\left( \omega \right) + m{k_{\textrm{aliasing}}} , \end{eqnarray}$$

(16)

$$\begin{eqnarray} k_n^{^{\prime\prime}}\left( \omega \right) = - {k_n}\left( \omega \right) + m{k_{\textrm{aliasing}}} , \end{eqnarray}$$

(17)

$$\begin{eqnarray} k_n^{^{\prime\prime\prime}}\left( \omega \right) = m{k_{\textrm{aliasing}}} . \end{eqnarray}$$

(18)

Where |$m = 0, \pm 1, \cdots $|⁠, |${k_{\textrm{aliasing}}} = 2\pi /\Delta x$| is the aliasing wavenumber. Their corresponding expressions for phase velocity are

$$\begin{eqnarray} c_n^{\prime}\left( \omega \right) = \frac{\omega }{{k_n^{\prime}\left( \omega \right)}} = \frac{{{c_n}\left( \omega \right)}}{{1 + m{k_{\rm{aliasing}}}{c_n}\left( \omega \right)/\omega }}\ , \end{eqnarray}$$

(19)

$$\begin{eqnarray} c_n^{^{\prime\prime}}\left( \omega \right) = \frac{\omega }{{k_n^{^{\prime\prime}}\left( \omega \right)}} = \frac{{{c_n}\left( \omega \right)}}{{ - 1 + m{k_{\rm{aliasing}}}{c_n}\left( \omega \right)/\omega }}\ , \end{eqnarray}$$

(20)

$$\begin{eqnarray} c_n^{^{\prime\prime\prime}}\left( \omega \right) = \frac{\omega }{{k_n^{^{\prime\prime\prime}}\left( \omega \right)}} = \frac{\omega }{{m{k_{\textrm{aliasing}}}}} . \end{eqnarray}$$

(21)

Fig. 6(a) shows the actual wavenumber |${k_n}( \omega )$| (⁠|$m = 0$|⁠) and several alias |${k_n}( \omega ) + m{k_{\textrm{aliasing}}}$|with |$m = - 2, - 1,1,2$| in the f–k domain. Their corresponding pattern in the f–v domain are shown in Figs 6(b)–(e). The mode branches for |$m = 0$| in Fig. 6(d) is the actual mode. The higher-order aliasing with larger m, for example, the alias for |$m = 2$| and |$m = 1$| shown in Figs 6(b) and (c), have a low velocity. They usually have less impact on the identification of the actual mode branches since they are located in the bottom of the dispersion image (see Figs 7a and c for m = 1 and m = 2). The aliasing for |$m = - 1$| can be divided into two parts (Fig. 6e). The part with negative velocity tends to infinity. The other part located in the upper right region has a positive velocity and has the potential to be misidentified as the high mode branches (see Figs 7a and c for |$m = - 1$|⁠).

$The illustration of the different types of the aliasing. (a) The aliasing $k = {k_n} + m{k_{\textrm{aliasing}}}$ associated to the positive wavenumber in the f–k domain. The corresponding results in the f–v domain are presented in panel (b) for $m = 2$, (c) for $m = 1$, (d) for $m = 0$ and (e) for $m = - 1$. (f) The aliasing $k = - {k_n} + m{k_{\textrm{aliasing}}}$ associated to the negative wavenumber in the f–k domain. The corresponding results in the f–v domain are presented in (g) for $m = 2$, (h) for $m = 1$, (i) for $m = 0$ and (j) for $m = - 1$. The actual mode branches are the result shown in (d) for m = 0. The black dashed lines denote $m{k_{\textrm{aliasing}}}$ or $- m{k_{\textrm{aliasing}}},m = \pm 1, \pm 2$.$

Figure 6.

The illustration of the different types of the aliasing. (a) The aliasing |$k = {k_n} + m{k_{\textrm{aliasing}}}$| associated to the positive wavenumber in the f–k domain. The corresponding results in the f–v domain are presented in panel (b) for |$m = 2$|⁠, (c) for |$m = 1$|⁠, (d) for |$m = 0$| and (e) for |$m = - 1$|⁠. (f) The aliasing |$k = - {k_n} + m{k_{\textrm{aliasing}}}$| associated to the negative wavenumber in the f–k domain. The corresponding results in the f–v domain are presented in (g) for |$m = 2$|⁠, (h) for |$m = 1$|⁠, (i) for |$m = 0$| and (j) for |$m = - 1$|⁠. The actual mode branches are the result shown in (d) for m = 0. The black dashed lines denote |$m{k_{\textrm{aliasing}}}$| or |$- m{k_{\textrm{aliasing}}},m = \pm 1, \pm 2$|⁠.

Open in new tab Download slide

$The theoretical dispersion curves and typical artifacts for the layered model in Table 1. The black solid lines denote the theoretical dispersion curves. Panels (a) and (b) are respectively the aliasing associated to ${k_n} + m{k_{\textrm{aliasing}}}$ and $- {k_n} + m{k_{\textrm{aliasing}}}$ for the interval $\Delta x = 15$ m. Panels (c) and (d) are the same as (a) and (b) but for the interval $\Delta x = 15$ m. The dashed lines represent the aliasing associated with different m labelled with the corresponding colours in the legend.$

Figure 7.

The theoretical dispersion curves and typical artifacts for the layered model in Table 1. The black solid lines denote the theoretical dispersion curves. Panels (a) and (b) are respectively the aliasing associated to |${k_n} + m{k_{\textrm{aliasing}}}$| and |$- {k_n} + m{k_{\textrm{aliasing}}}$| for the interval |$\Delta x = 15$| m. Panels (c) and (d) are the same as (a) and (b) but for the interval |$\Delta x = 15$| m. The dashed lines represent the aliasing associated with different m labelled with the corresponding colours in the legend.

Open in new tab Download slide

Figs 6(f)–(j) show the aliasing associated with |$- {k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| in the f–k and f–v domain, respectively. The aliasing for |$m \le 0$| has the negative velocity (e.g. the aliasing for |$m = 0$| and |$m = - 1$| in Figs 6(i) and (j) and it does not affect the dispersion image. The aliasing of |$- {k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| at |$m = 1$| exhibit mirror symmetry with those of |${k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| at |$m = - 1$|⁠. The aliasing with positive velocity for |$m \ge 1$| crosses the actual mode branches and known as the crossed artifacts (see Figs 7b and d for m = 1 and m = 2).

The dispersion curves and aliasing we finally observed in a typical dispersion image are shown in Fig. 7, where the theoretical dispersion curves and the low-order artifacts are present. The black solid lines in Fig. 7 denote the theoretical dispersion curves. Figs 7(a) and (b) show the aliasing associated with |${k_n} + m{k_{\textrm{aliasing}}}$| and |$- {k_n} + m{k_{\textrm{aliasing}}}$| (⁠|$m = - 1,1,2$|⁠) for |$\Delta x = 15$| m. The results for|$\Delta x = 30$| m are presented in Figs 7(c) and (d). The radial artifacts associated with |$m{k_{\textrm{aliasing}}}$| are indicated by the blue solid lines in Figs 7(b) and (d).

The aliasing associated to |${k_n} + m{k_{\textrm{aliasing}}}$| with |$m > 0$| occurs at the bottom with low velocity and does not significantly affect the recognition of mode branches, as shown in Fig. 7(a) for m = 1 and m = 2. However, their location depends on the station spacing. For a larger station spacing, for example, |$\Delta x = 30$|⁠, the aliasing associated with |$m = 1$| crosses with the theoretical dispersion curves, as shown in Fig. 7(c). Similarly, the aliasing associated with |$m = - 1$| also depends on the station spacing. In Fig. 7(a), this aliasing is mainly located in the upper right region and far from the mode branches of interest. But in Fig. 7(c), it moves towards the lower frequency for |$\Delta x = 30$| and is likely to be recognized as the actual mode branch due to their similar pattern. The feature that the velocity of the aliasing tends to infinity provides a way of distinguishing them from the actual mode. In addition, once we have identified the actual mode branches, the position of the artifacts for |$m = - 1$| can be estimated using eq. (19), which offers a way to distinguish and reduce them from the actual mode. Such aliasing can also be reduced by random sampling on the wavefield, as discussed in section 6.

As for the aliasing associated with |$- {k_n} + m{k_{\textrm{aliasing}}}$| and |$m{k_{\textrm{aliasing}}}$| shown in Figs 7(b) and (d), they can be removed by changing the imaging conditions and simple data pre-processing. These aliasing are often referred to as ‘crossed’ and ‘radial’ artifacts (e.g. Luo et al. 2022; Cheng et al. 2023). The location and the number of crossed artifacts depends on the array configuration. It likely always crosses with the actual modes. The radial artifacts associated with |$m{k_{\textrm{aliasing}}}$| are strongly controlled by the minimum station spacing. If the station spacing is small, it appears mainly in the region around zero velocity at the bottom of the dispersion image, as shown in Fig. 7(b) for |$\Delta x = 15$|⁠. For a larger station interval, as shown in Fig. 7(d) for |$\Delta x = 30$|⁠, it will go through the actual mode branch, affecting the clarity and resolution of the dispersion image.

3.4 Removal of crossed artifacts based on Kramers–Kronig relations

According to the SI theory, the NCF retrieved from the seismic noise in the time domain is equivalent to the spatial autocorrelation coefficient in the frequency domain (e.g. Yokoi & Margaryan 2008; Tsai & Moschetti 2010; Lu 2021), which can be expressed as |${J_0}(kr)$|⁠. The spatial wavefield described by |${J_0}(kr)$| can be thought of as a standing wave with cylindrical wavefront, which is the superposition of two cylindrical travelling waves propagating inwards and outwards. The |${J_0}({k_n}r)$| can be expressed as:

$$\begin{eqnarray} {J_0}\left( {{k_n}r} \right) = 1/2\left[ {H_0^{\left( 1 \right)}\left( {{k_n}r} \right) + H_0^{\left( 2 \right)}\left( {{k_n}r} \right)} \right] , \end{eqnarray}$$

(22)

where |$H_0^{(1)}$| and |$H_0^{(2)}$| are respectively the first and second kind of Hankel function with zero-order. According to the convention of Fourier transform given in Appendix A, |$H_0^{(1)}$|represents the waves propagating inwards while |$H_0^{(2)}$| represents the waves propagating outwards. For the imaging condition in eq. (4), a plane wave represented by a cosine function is used to fit the standing wave via delay-and-sum. This means that two cylindrical waves propagating in opposite directions are fitted simultaneously. The wavenumbers could be also estimated for the waves propagating inwards, that is, propagating along direction |$- {k_n}$|⁠. The aliasing wavenumber |$k_n^{\prime}$| can be observed at |$k_n^{\prime} = - {k_n} + {k_{\textrm{aliasing}}}$| due to the symmetry with the target wavenumber at |$+ {k_n}$| direction.

Therefore, if only the wavefield propagating in |$+ {k_n}$| direction is considered in the imaging condition, the aliasing associated with |$- {k_n} + {k_{\textrm{aliasing}}}$| will be eliminated. Forbriger (2003) applied this idea to process the seismogram gathered from an active source, where the Bessel function is replaced by the Hankel function in the wavefield transformation. This idea is extended to F-J method of the virtual records represented by NCFs in which the Hilbert transform is used to construct the complete GF and obtain the propagating waves in one direction (Xi et al. 2021; Zhou & Chen, 2022). In this paper, we present a new BF imaging condition for removal of crossed artifacts. Different from the modified F-J method, the construction of the complete GF is not necessary. All the variables involved in the proposed new BF imaging conditions are the elements in CSDM.

4 NEW IMAGING CONDITIONS OF BEAMFORMING

4.1 Causal and acausal parts of GF

The theoretical basis of the new proposed imaging conditions is the SI theory. In this section, we briefly revisit the SI theory based on the plane wave model (e.g. Boschi et al. 2013; Lu 2021).

As shown in the top panel of Fig. 8(a), a monochromatic plane wave is incident onto stations A and B from the direction |$\theta $|⁠. |${{\bf{x}}_\mathrm{A}}$| and |${{\bf{x}}_\mathrm{B}}$| denote the location of A and B, respectively. The spatial coherence between stations A and B in the frequency domain can be expressed as |$\exp [ - \iota {{\bf{k}}_n} \cdot ({{\bf{x}}_\mathrm{B}} - {{\bf{x}}_\mathrm{A}})]$|⁠. Their real and imaginary parts are shown in the panels of the second and third row of Fig. 8(a). Fourier integration of the results for monochromatic plane waves over a certain bandwidth gives results for the transient case, as shown schematically in the bottom panels. Because the source is not located on the line connecting two stations, the measured velocity is the apparent one and the lagged time of the waveforms (solid line) is not the position (dashed line) that the GF would have appeared.

$Schematic diagram for the reconstruction of Green's function by seismic interferometry. The panels in the top row show the distribution of the noise sources radiating a far-field monochromatic plane wave. The panels in the middle two rows show the real and imaginary parts of the spatial coherency of a monochromatic plane wave. The panels in the bottom row show the sketch of the cross-correlation function between two stations in the time domain. The left column (a) shows the case of a single source. The middle (b) and right (c) columns show the case where the sources are located in the azimuth range of $\theta \in ( {\pi ,2\pi } )$ and $\theta \in ( {0,\pi } )$, respectively. $\theta $ represent the azimuth of the source, as shown by the top panel in (a). The grey area denotes the Fresnel region. ${r_{\mathrm{ AB}}}$ represents the distance between station ${{\bf{x}}_\mathrm{A}}$ and station ${{\bf{x}}_\mathrm{B}}$.$

Figure 8.

Schematic diagram for the reconstruction of Green's function by seismic interferometry. The panels in the top row show the distribution of the noise sources radiating a far-field monochromatic plane wave. The panels in the middle two rows show the real and imaginary parts of the spatial coherency of a monochromatic plane wave. The panels in the bottom row show the sketch of the cross-correlation function between two stations in the time domain. The left column (a) shows the case of a single source. The middle (b) and right (c) columns show the case where the sources are located in the azimuth range of |$\theta \in ( {\pi ,2\pi } )$| and |$\theta \in ( {0,\pi } )$|⁠, respectively. |$\theta $| represent the azimuth of the source, as shown by the top panel in (a). The grey area denotes the Fresnel region. |${r_{\mathrm{ AB}}}$| represents the distance between station |${{\bf{x}}_\mathrm{A}}$| and station |${{\bf{x}}_\mathrm{B}}$|⁠.

Open in new tab Download slide

As shown in the top panel of Fig. 8(b), if the sources are continuously distributed in |$( {\pi ,2\pi } )$|with unit amplitude, |${C_{\mathrm{ AB}}}( \omega )$| between stations A and B can be calculated by integrating |$\exp [ - \iota {{\bf{k}}_n} \cdot ({{\bf{x}}_\mathrm{B}} - {{\bf{x}}_\mathrm{A}})]$| for those single sources. It can be expressed as (e.g. Lu 2021)

$$\begin{eqnarray} {C_{\mathrm{ AB}}}\left( \omega \right) = {J_0}\left( {k{r_{\mathrm{ AB}}}} \right) - \iota {S_0}\left( {k{r_{\mathrm{ AB}}}} \right) \end{eqnarray}$$

(23)

where |${S_0}$| is the zero-order Struve function. k is the scalar wavenumber and |${r_{\mathrm{ AB}}}$| is the distance between stations A and B. The real and imaginary parts of eq. (23) are shown in the panels of the second and third row of Fig. 8(b). The corresponding transient solutions are schematically shown in the bottom panel, representing the waves propagating from station A to B.

Similarly, for the case shown in the top panel of Fig. 8(c), where the waves incident from|$\theta \in ( {0,\pi } )$|⁠, |${C_{\mathrm{ AB}}}( \omega )$| can be expressed as (e.g. Lu 2021)

$$\begin{eqnarray} {C_{\mathrm{ AB}}}\left( \omega \right) = {J_0}\left( {k{r_{\mathrm{ AB}}}} \right) + \iota {S_0}\left( {k{r_{\mathrm{ AB}}}} \right) \end{eqnarray}$$

(24)

The bottom three panels of Fig. 8(c) shows the results in the frequency domain and the corresponding transient solution. It represents the waves propagating from station B to A. If the sources are uniformly distributed in all azimuth, |${C_{\mathrm{ AB}}}( \omega )$| is the sum of eqs (23) and (24), i.e. the well-known result |${J_0}(k{r_{\mathrm{ AB}}})$|⁠, denoting a standing wave with a cylindrical wavefront.

According to eqs (23) and (24), we have already obtained the propagating wave in one direction, and what we need to do is directly utilize the causal or acausal parts of NCFs to obtain the propagating wave in one direction. The imaginary part of the causal NCFs corresponds to |${S_0}( {kr} )$|⁠. In the far field, it tends to zero-order Bessel function |${Y_0}( {kr} )$| of the second kind. In the near field (⁠|$r \to 0$|⁠), |${Y_0}( {kr} )$| is infinite but |${S_0}( {kr} )$| is finite (see eqs 9.1.13, 12.1.4, 12.1.30 and figs 9.1, 12.1 in Abramowitz & Stegun 1964).

To consider the propagating waves in one direction, we construct the CSDM using the causal part of the NCFs by

$$\begin{eqnarray} C_{ij}^{\textrm{causal}}\left( \omega \right) = {\rm{FT}}\left\{ {\left[ {\begin{array}{@{}*{4}{c}@{}} 0&{\mathrm{ NCF}({r_{12}},t)}& \cdots &{\mathrm{ NCF}({r_{1N}},t)}\\ {\mathrm{ NCF}({r_{21}}, - t)}&0&{}& \vdots \\ \vdots &{}&0&{\mathrm{ NCF}({r_{(N - 1)N}},t)}\\ {\mathrm{ NCF}({r_{N1}}, - t)}& \cdots &{\mathrm{ NCF}({r_{N(N - 1)}}, - t)}&0 \end{array}} \right]} \right\} \\ \end{eqnarray}$$

(25)

where N represents the number of stations in the array. The subscripts i and j denote the indices of the stations. |${\rm{FT}}\{ {} \}$| represents the Fourier transform.|$\mathrm{ NCF}({r_{ij}},t)$| and |$\mathrm{ NCF}({r_{ij}}, - t)$| represent the time-domain field propagating from station i to station j. |$\mathrm{ NCF}({r_{ij}},t)$| is the causal part, while |$\mathrm{ NCF}({r_{ij}}, - t)$| is the acausal part. According to the reciprocity theorem, |$\mathrm{ NCF}({r_{ij}}, - t)$| would represent the causal part (⁠|$\mathrm{ NCF}({r_{ji}},t)$|⁠) propagating from station j to station i. Thereby, both |$\mathrm{ NCF}({r_{ij}},t)$| and |$\mathrm{ NCF}({r_{ji}}, - t)$| represent the waves propagating from station i to station j. Only waves propagating in one direction is considered in |$C_{ij}^{\textrm{causal}}( \omega )$|⁠. Based on |$C_{ij}^{\textrm{causal}}( \omega )$| rather than the original |${C_{ij}}( \omega )$|⁠, the new imaging condition can be constructed to remove the crossed artifacts.

4.2 New imaging conditions of MCBF for removal of crossed artifacts

The new imaging condition of MCBF shown in eq. (4) is expressed as:

$$\begin{eqnarray} {\overline {\textrm{MCBF}} _{\textrm{imag}}}\left( {k,\omega } \right) = \frac{1}{2}\left( {\mathrm{ CC} + \mathrm{ SS}} \right) \end{eqnarray}$$

(26)

with

$$\begin{eqnarray} \mathrm{ CC} = \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {\sqrt {\pi {r_{ij}}k} {\mathop{\rm{Re}}\nolimits} \left[ {C_{ij}^{\textrm{causal}}\left( \omega \right)} \right]\cos \left( {k{r_{ij}} - \frac{\pi }{4}} \right)} } , \end{eqnarray}$$

(27)

$$\begin{eqnarray} \mathrm{ SS} = \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {\sqrt {\pi {r_{ij}}k} {\mathop{\rm{Im}}\nolimits} \left[ { - C_{ij}^{\textrm{causal}}\left( \omega \right)} \right]\sin \left( {k{r_{ij}} - \frac{\pi }{4}} \right)} } . \end{eqnarray}$$

(28)

Eq. (27) is equivalent to the original imaging condition shown in eq. (4). In the far field, the real and imaginary part of eq. (23) can be approximated by |$\sqrt {{2 {/ {\vphantom {2 {\pi kr}}} } {\pi kr}}} \cos (kr - {\pi {/ {\vphantom {\pi 4}} } 4})$| and |$\sqrt {{2 {/ {\vphantom {2 {\pi kr}}} } {\pi kr}}} \sin (kr - {\pi {/ {\vphantom {\pi 4}} } 4})$| (eqs 9.2.1 and 12.1.34 in Abramowitz & Stegun 1964), respectively, and eqs (27) and (28) can be written as (see Appendix B for details)

$$\begin{eqnarray} \mathrm{ CC} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left[ {\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} + \sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}} \right] \\ \mathrm{ SS} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left[ {\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} - \sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}} \right] \\ \end{eqnarray}$$

(29)

where |${\rm{sinc(}}x) = {{\sin x} {/ {\vphantom {{\sin x} x}} } x}$| is the sinc function, and |$M = {{N(N - 1)} {/ {\vphantom {{N(N - 1)} 2}} } 2}$|⁠. R the maximum interstation distance of the array. Eq. (29) shows the terms containing only the argument |${k_n} - k$| are separated from the terms containing only argument |${k_n} + k$|⁠. Eq. (26) can be analytically written as

$$\begin{eqnarray} {\overline {\textrm{MCBF}} _{\textrm{imag}}}\left( {k,\omega } \right) = \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} . \end{eqnarray}$$

(30)

The crossed artifacts can be estimated using

$$\begin{eqnarray} \mathrm{ MCB}{\mathrm{ F}_{\textrm{alias}}}\left( {k,\omega } \right) = \frac{1}{2}\left( {\mathrm{ CC} - \mathrm{ SS}} \right) \end{eqnarray}$$

(31)

and which can be analytically written as (see Appendix B for details)

$$\begin{eqnarray} \mathrm{ MCB}{\mathrm{ F}_{\textrm{alias}}}\left( {k,\omega } \right) = \frac{{2\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}. \\ \end{eqnarray}$$

(32)

A simple numerical example is used to illustrate how the new imaging conditions eliminate the crossed artifacts. It is assumed that the cylindrical wavefield described by |${J_0}(kr)$|contains three dimensionless wavenumbers |${k_1} = 1$|⁠, |${k_2} = 10$| and |${k_3} = 25$| with unit amplitude. The |$C_{ij}^{\textrm{causal}}( \omega )$| for such a wavefield can be expressed as

$$\begin{eqnarray} C_{ij}^{\textrm{causal}}\left( \omega \right) = {J_0}({k_1}{r_{ij}}) + {J_0}({k_2}{r_{ij}}) + {J_0}({k_3}{r_{ij}}) \end{eqnarray}$$

(33)

In simulation, the spatial wavefield is uniformly sampled with |$\Delta r = 0.2$| and the maximum distance |${r_{\max }} = 8.0$|⁠. We investigate the effect of different imaging conditions on the result by comparing the beamformer output of such an ideal cylindrical wavefield by substituting eq. (33) into eqs (27) and (28).

Fig. 9(a) shows the beamformer output of eq. (27), the original imaging condition. Two aliasing wavenumbers (⁠|$k_1^{\prime} = - {k_1} + 2{k_{{\rm{aliasing}}}} = 30.4$|⁠, |$k_2^{\prime} = - {k_2} + 2{k_{{\rm{aliasing}}}} = 21.4$|⁠) caused by the symmetry of |${k_1}$| and |${k_2}$| are observed at |$k > {k_{{\rm{aliasing}}}}$|⁠. The aliasing wavenumber associated with |${k_3}$| appears at |$k_3^{\prime} = - {k_3} + {k_{{\rm{aliasing}}}} = 6.4$| due to the symmetry and finite sampling. The target wavenumbers appear at the positive peaks of the beamformer output, shaped as |$\cos ( {k - {k_n}} ){\rm{sinc}}( {k - {k_n}} )$|⁠. The aliasing wavenumbers |$k_n^{\prime}$| (n = 1, 2, 3) appear at the zero-crossing points between two positive and negative peaks, locally shaped as |$\sin ( {k + {k_n}} ){\rm{sinc}}( {k + {k_n}} )$|⁠.

$An illustration for aliasing wavenumbers and its reduction. The results are calculated by different imaging conditions for a cylindrical wavefield containing three horizontal wavenumbers ${k_1} = 1$, ${k_2} = 10$ and ${k_3} = 25$, which are indicated by vertical orange dashed lines. The aliasing wavenumbers $k_1^{\prime}$, $k_2^{\prime}$ and $k_3^{\prime}$ are denoted by vertical green dashed lines. ${k_{{\rm{aliasing}}}} = \pi /{r_{\min }} = 15.7$ is denoted by grey vertical dashed line. (a) The result of eq. (4), that is, the CC of eq. (27). (b) The result of SS of eq. (28). (c) The result of the new imaging condition shown in eq. (26). (d)The aliasing wavenumber obtained by eq. (31).$

Figure 9.

An illustration for aliasing wavenumbers and its reduction. The results are calculated by different imaging conditions for a cylindrical wavefield containing three horizontal wavenumbers |${k_1} = 1$|⁠, |${k_2} = 10$| and |${k_3} = 25$|⁠, which are indicated by vertical orange dashed lines. The aliasing wavenumbers |$k_1^{\prime}$|⁠, |$k_2^{\prime}$| and |$k_3^{\prime}$| are denoted by vertical green dashed lines. |${k_{{\rm{aliasing}}}} = \pi /{r_{\min }} = 15.7$| is denoted by grey vertical dashed line. (a) The result of eq. (4), that is, the CC of eq. (27). (b) The result of SS of eq. (28). (c) The result of the new imaging condition shown in eq. (26). (d)The aliasing wavenumber obtained by eq. (31).

Open in new tab Download slide

Fig. 9(b) shows the result of SS in eq. (28). The target wavenumbers also appear at the positive peaks. The aliasing wavenumbers appear at the zero-crossing points but it has the opposite phase with CC generated by eq. (27). Therefore, the aliasing wavenumbers disappear by summing them under the new imaging condition, as shown in Fig. 9(c), where only the target wavenumbers |$k = {k_n}$| propagating in |${\bf{k}}$| direction remain. Similarly, the result of subtracting SS from CC gives the aliasing wavenumbers, as shown in Fig. 9(d), where only the aliasing wavenumbers due to the wave propagating in |$- {\bf{k}}$| direction remain.

4.3 New imaging conditions of WCBF

For WCBF, the new imaging conditions can be expressed as

$$\begin{eqnarray} \begin{array}{@{}l@{}} \mathrm{ WCB}{\mathrm{ F}_{\textrm{imag}}}\left( {k,\omega ,\theta } \right) = \left( {\mathrm{ CC} + \mathrm{ SS}} \right)/2\\ \mathrm{ WCB}{\mathrm{ F}_{\textrm{alias}}}\left( {k,\omega ,\theta } \right) = \left( {\mathrm{ CC} - \mathrm{ SS}} \right)/2 \end{array} \end{eqnarray}$$

(34)

where

$$\begin{eqnarray} \begin{array}{@{}l@{}} \mathrm{ CC} = \frac{1}{{N{{\left( \theta \right)}^2}}}\sum\limits_i^{N\left( \theta \right)} {\sum\limits_j^{N\left( \theta \right)} {{\mathop{\rm{Re}}\nolimits} \left[ {C_{ij}^{\textrm{cause}}\left( {\omega ,{\theta _{ij}}} \right)} \right]\cos \left[ {k{r_{ij}}\cos \left( {{\theta _{ij}} - \theta } \right)} \right]k{r_{ij}}} } \\ \mathrm{ SS} = \frac{1}{{N{{\left( \theta \right)}^2}}}\sum\limits_i^{N\left( \theta \right)} {\sum\limits_j^{N\left( \theta \right)} { - {\mathop{\rm{Im}}\nolimits} \left[ {C_{ij}^{\textrm{cause}}\left( {\omega ,{\theta _{ij}}} \right)} \right]\sin \left[ {k{r_{ij}}\cos \left( {{\theta _{ij}} - \theta } \right)} \right]k{r_{ij}}} } \\ \ ,\theta - \pi /2 \le {\theta _{ij}} \le \theta + \pi /2 \end{array} \end{eqnarray}$$

(35)

The projection operation is performed only for |${\theta _{ij}} \in [ {\theta - \pi /2,\theta + \pi /2} ]$|⁠. |$\mathrm{ WCB}{\mathrm{ F}_{\textrm{imag}}}( {k,\omega ,\theta } )$| gives the true wavenumber and |$\mathrm{ WCB}{\mathrm{ F}_{\textrm{alias}}}( {k,\omega ,\theta } )$| gives the aliasing associated with negative wavenumbers.

5 EXAMPLES AND DISCUSSION

5.1 Synthetic data

In this section, we investigate the extraction of multimode dispersion curves of Rayleigh wave using new imaging conditions based on the synthetic data. Fig. 10(a) shows the configuration of the array and sources used in the simulation. The 10 000 point sources randomly distributed over a 1.0–1.5 km annular region. The receiver array consists of 144 stations, which are distributed uniformly within an area of 0.2 km × 0.2 km (Fig. 10b). The station spacing |$\Delta d$| is 0.015 km. Fig. 10(c) shows the ARF. The square array with equal spacing produces an ARF with infinite subpeaks that have the same magnitude as the main one. Fig. 10(d) shows the cross-section of ARF every two degrees. The subpeaks give the wavenumbers that generate the aliasing, as already seen in Section 3.2. The aliasing wavenumber associated with the first and second subpeak are |$k_{\textrm{aliasing}}^1 = 2\pi /\Delta d = 419$| and |$k_{\textrm{aliasing}}^2 = 2\pi /\Delta d \cdot \sqrt 2 = 592$|⁠, respectively.

Figure 10.

(a) The configuration of the receiver array (blue) and sources (orange). (b) The distribution of the receiver array. (c) The ARF of the receiver array. (d) The cross-sections of ARF in radial direction every two degrees. The subpeaks and the corresponding wavenumber value are labled in (d).

Open in new tab Download slide

The layered model shown in Table 1 is considered again. The method to synthesize the data is the same as Qin & Lu (2024). The continuous noise records are synthesized by calculating the vertical velocity component of the Rayleigh wave GF, assuming that the intensity of the point source is random and the central frequencies of the source functions (Ricker wavelet) distributed randomly at 0–25 Hz. The NCFs are then computed by cross-correlating the synthetic noise recordings (e.g. Bensen et al. 2007). Fig. 11(a) shows the NCFs for some station-pairs. With increasing distance, the wave packet associated with higher modes are gradually separated. As shown in Fig. 11(b), the upper triangular matrix of CSDM is constructed using the Fourier transform of the causal part of the NCFs. The lower one is formed by the Fourier transform of the time-reversed version of the acausal part of NCFs.

Figure 11.

An illustration to construct the CSDM using the waves propagating in one direction. (a) The synthetic NCFs, filtered in 5–20 Hz. (b) Schematic to construct the CSDM using NCFs.

Open in new tab Download slide

The elements in CSDM constitute a data set of the sampling on the spatial wavefield. Figs 12(a) and (b) show the real and imaginary parts of the wavefield at 5 and 7.5 Hz, respectively. It shows the regular array only samples the wavefield at specific positions. For the laterally isotropic structure, the wavefield at different azimuths can be projected onto the radial direction, as shown in the bottom panels of Fig. 12. Radially, the real and imaginary parts of the wavefield can be represented by |$\sum\limits_n {{A_n}{J_0}({k_n}r)} $| and |$\sum\limits_n {{A_n}{S_0}({k_n}r)} $|⁠, respectively, where |${A_n}$| denotes the amplitude of mode n. At 5 Hz, only the fundamental mode is generated. The real and imaginary parts fit well with |${J_0}({k_0}r)$| and |${S_0}({k_0}r)$|⁠, respectively. At 7.5 Hz, both the fundamental mode and the first overtone are generated. The real and imaginary parts approximately fit with |${A_0}{J_0}({k_0}r) + {A_1}{J_0}({k_1}r)$| and |${A_0}{S_0}({k_0}r) + {A_1}{S_0}({k_1}r)$|⁠, respectively.

$The synthetic wavefield sampled by the regular array (see Appendix C for sampling positions). (a) The real part of the wavefield at 5 and 7.5 Hz. The 2-D spatial distribution is shown in the upper panel. The radially distributed wavefield is presented in the lower panel. The blue lines represent ${J_0}( {{k_0}r} )$ for 5 Hz and ${A_0}{J_0}( {{k_0}r} ) + {A_1}{J_0}( {{k_1}r} )$ for 7.5 Hz. (b) The imaginary part of the wavefield. The orange lines represent${S_0}( {{k_0}r} )$ for 5 Hz and ${A_0}{S_0}( {{k_0}r} ) + {A_1}{S_0}( {{k_1}r} )$ for 7.5 Hz.$

Figure 12.

The synthetic wavefield sampled by the regular array (see Appendix C for sampling positions). (a) The real part of the wavefield at 5 and 7.5 Hz. The 2-D spatial distribution is shown in the upper panel. The radially distributed wavefield is presented in the lower panel. The blue lines represent |${J_0}( {{k_0}r} )$| for 5 Hz and |${A_0}{J_0}( {{k_0}r} ) + {A_1}{J_0}( {{k_1}r} )$| for 7.5 Hz. (b) The imaginary part of the wavefield. The orange lines represent|${S_0}( {{k_0}r} )$| for 5 Hz and |${A_0}{S_0}( {{k_0}r} ) + {A_1}{S_0}( {{k_1}r} )$| for 7.5 Hz.

Open in new tab Download slide

Fig. 13 shows the result of CC, SS, (CC + SS)/2 and (CC − SS)/2 in the f–k and f–v domain. In terms of the generation of the azimuthally averaged dispersion image, the results obtained by WCBF and MCBF are equivalent (Qin & Lu 2024). Only the results of MCBF are shown. As expected, the crossed artifacts are removed using the new imaging condition, as shown by the result of (CC + SS)/2. As indicated by black arrows, the aliasing associated with |${k_n} + m{k_{\textrm{aliasing}}}$| for m = 1 is not removed, but it can be predicted by eq. (19), as discussed in Section 3.3.

$The dispersion image generated by different imaging conditions for the synthetic data of a regular array. The solid black lines denote the locations of the wavenumber $k_{\textrm{aliasing}}^1$, $k_{\textrm{aliasing}}^2$ and $k_{\textrm{aliasing}}^3$. The black dashed lines denote the location of $k_{\textrm{aliasing}}^1/2$ and $k_{\textrm{aliasing}}^2/2$. The crossed artifacts appear in CC are removed by the new imaging condition of (CC + SS)/2. The black arrows denote the aliasing associated with ${k_n} + m{k_{\textrm{aliasing}}}$, which cannot be removed by new imaging condition but can be predicted by eq. (19).$

Figure 13.

The dispersion image generated by different imaging conditions for the synthetic data of a regular array. The solid black lines denote the locations of the wavenumber |$k_{\textrm{aliasing}}^1$|⁠, |$k_{\textrm{aliasing}}^2$| and |$k_{\textrm{aliasing}}^3$|⁠. The black dashed lines denote the location of |$k_{\textrm{aliasing}}^1/2$| and |$k_{\textrm{aliasing}}^2/2$|⁠. The crossed artifacts appear in CC are removed by the new imaging condition of (CC + SS)/2. The black arrows denote the aliasing associated with |${k_n} + m{k_{\textrm{aliasing}}}$|⁠, which cannot be removed by new imaging condition but can be predicted by eq. (19).

Open in new tab Download slide

5.2 Field data

The results generated by the new imaging condition for the array shown in Fig. 1(a) is given in Fig. 14. In comparison with Fig. 1(c), the crossed artifacts associated with |$- {k_n} + m{k_{\textrm{aliasing}}}$| are removed in Fig. 14(b). The artifacts associated with |${k_n} + m{k_{\textrm{aliasing}}}$| are characterized similarly to the actual mode branch. If their amplitude is comparable to that of the actual mode, they would affect the mode recognition. Fortunately, once we get the phase velocity of the two actual modes, we can predict the location of these artifacts based on eq. (19). For the array shown in Fig. 1(a), two minimum aliasing wavenumbers are |$k_{\textrm{aliasing}}^1 = 6.3$| and |$k_{\textrm{aliasing}}^2 = 8.9$|⁠. The predicted artifacts for m = 1, m = −1, and m = −2 are presented in Fig. 14(c). The estimation of their locations would help to identify them from the actual modes.

$The beamformer output and dispersion image generated by new imaging condition for the array in Fig. 1(a). (a) Beamformer output at 0.4, 0.5 and 0.7 Hz. (b) Azimuth-averaged dispersion image generated by new imaging condition. (c) The mode branches of the picked fundamental mode (black solid line) and the first overtone. The coloured solid lines denote the aliasing ${k_n} + m{k_{\textrm{aliasing}}}$ for m = 1, m = −1, and m = −2 associated with the fundamental mode, while the coloured dashed lines denote those associated with the first overtone. These predicted artifacts are also presented in (b) with black lines. (d) The comparison of the picked mode branches from the dispersion image generated by original and new imaging condition. Some slight fluctuations are marked by the black arrows.$

Figure 14.

The beamformer output and dispersion image generated by new imaging condition for the array in Fig. 1(a). (a) Beamformer output at 0.4, 0.5 and 0.7 Hz. (b) Azimuth-averaged dispersion image generated by new imaging condition. (c) The mode branches of the picked fundamental mode (black solid line) and the first overtone. The coloured solid lines denote the aliasing |${k_n} + m{k_{\textrm{aliasing}}}$| for m = 1, m = −1, and m = −2 associated with the fundamental mode, while the coloured dashed lines denote those associated with the first overtone. These predicted artifacts are also presented in (b) with black lines. (d) The comparison of the picked mode branches from the dispersion image generated by original and new imaging condition. Some slight fluctuations are marked by the black arrows.

Open in new tab Download slide

Compared to Fig. 1(b), the dispersion image generated by the new imaging condition as a higher resolution due to the removal of the crossed artifacts. As a result, with the same measurement uncertainty, phase velocities can be extracted over a wider bandwidth using new imaging condition, as shown in Fig. 14(c). Moreover, the phase velocities extracted from Fig. 14(b) are smoother along the mode branches, and the slight fluctuations in the velocities extracted from Fig. 1(b) due to crossed artifacts disappear.

6 RANDOM SAMPLING AND ARRAY DESIGN

The aliases arise from the regular sampling on the spatial wavefield (e.g. Zwartjes & Sacchi 2007; Snider & Wakin 2022). The crossed artifacts associated with |$- {k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| can be removed using the new imaging condition. However, as indicated by the black arrows in Fig. 13, the aliases associated with |${k_n} + m{k_{\textrm{aliasing}}}$| remain. In this section, we show these aliases can be reduced by replacing the regular array with the random one.

6.1 Random sampling helps reduce aliasing: synthetic and field data

Fig. 15(a) shows the configuration of a random array consisting of 144 stations. The number of stations and the aperture are the same as in Fig. 10(b). For the given number of stations and the covering area, we use the pseudo-random number generator (PCG-64, O'Neill 2014) to generate the positions of the station. Fig. 15(b) shows the ARF of the array. Fig. 15(c) shows the cross-sections of ARF every two degrees. In comparison with Figs 10(c) and (d), the amplitude of the subpeaks is much smaller than the main one. Less aliasing is therefore expected in the dispersion image for such a random array.

Figure 15.

(a) The station distribution of the random array. (b) The ARF of the array. (c) The cross-sections (grey lines) of ARF in radial direction every two degrees. The black dashed line presents the azimuth-averaged result.

Open in new tab Download slide

We perform the synthetic test for the random array and compare them with that for the regular one. The model and parameters used to synthesize the data are the same as those used for the regular array. Fig. 16 shows the spatial wavefield sampled by the random array at 5 and 7.5 Hz. It shows the wavefield is sampled adequately along the radial direction.

Figure 16.

The same as Fig. 12 but the wavefield is sampled by the random array shown in Fig. 15(a).

Open in new tab Download slide

Figs 17(a) and (b) show the dispersion images in the f–k and f–v domain generated by different imaging conditions. The aliases are significantly reduced compared to the result for the array with the regular configuration. Even for the result obtained by the original imaging condition, the artifacts are removed due to the random sampling. Moreover, the aliases indicated by the arrows in Fig. 13 have been eliminated in Fig. 17. This implies that random sampling can reduce all the aliases associated with the array configuration.

Figure 17.

Azimuth-averaged dispersion image generated by different imaging conditions for the random array. (a) The dispersion image in f–k domain, (b) The dispersion image in f–v domain. From left to right, it shows the result of CC, SS, (CC + SS)/2 and (CC − SS)/2., respectively.

Open in new tab Download slide

For the periodic sampling of the regular array, the aliases associated with the sampling interval is coherent. Random sampling, on the other hand, breaks the constructive interference induced by periodic sampling and turns the coherent aliasing into Gaussian noise (e.g. Herrmann 2010; Snider & Wakin 2022). As a result, the aliases associated with periodic sampling, including aliases associated with |$- {k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| and |${k_n}( \omega ) + m{k_{\textrm{aliasing}}}$|⁠, are reduced. The dispersion image generated by the wavefield sampled randomly has higher clarity and fewer artifacts.

This also explains that the artifacts in the dispersion image in Fig. 1(c) are insignificant since the real array is not perfect regular. To demonstrate the case for an array with more randomly distributed stations, we resampled 56 stations from Fig. 1(a) randomly, as denoted by the blue triangles in Fig. 18(a). These stations are complementary to those considered in Figs 1(a). Fig. 18(b) shows the ARF of this array. It shows the amplitude of the sub-peaks decays rapidly with increasing wavenumber. Compared with Fig. 2(b), the amplitude of the main lobe in Fig. 18(b) is significantly larger than those of the surrounding subpeaks. The dispersion image generated by eqs (26) and (27) are presented in Figs 18(c) and (d), respectively. The new imaging condition removes the crossed artifacts successfully. Furthermore, in contrast to Fig. 1(c), the crossed artifacts, as well as the other artifacts, are much weaker for the random array even for the original imaging conditions. This means random sampling plays a quite good role in reducing artifacts.

Figure 18.

An illustration of the dispersion image generated by original and new imaging conditions for a random array. (a) The station configuration of the array (blue triangles) (The stations denoted by the grey triangles are those considered in Fig. 1a). (b) ARF of the array. (c) The dispersion image generated by original imaging condition. (d) The dispersion image generated by new imaging condition. (e) The comparison of the picked mode branches from the dispersion image generated by original and new imaging condition.

Open in new tab Download slide

6.2 Random array design and performance comparison

Since aliasing is associated with the periodic sampling, it is natural to think of designing a seismic array that samples the wavefield non-periodically to reduce aliasing. Recently, an aperiodic monolith called ‘the Hat’ has been discovered that allows for aperiodic tiling in the plane (Smith et al., 2024). Mordret & Grushin (2024) discussed the possibility of introducing aperiodic monolith into the design of the seismic array. Assuming that the station is placed at the vertex of ‘the Hat’ family of tiles, they analyzed the ARF and found that the aperiodic monolith seismic array outperforms regular arrays.

In this section, we discuss the performance of several seismic arrays in generating high-resolution multimode dispersion image. These arrays are composed of the tiles of the Hat family or are designed using random sampling schemes. A feasible array design scheme is then suggested.

Figs 19(a)–(c) show three typical shapes of ‘the Hat’ family: the Hat, the Specter, and the Turtle (Mordret & Grushin 2024; Smith et al. 2024), which have been shown to be capable of tiling the plane aperiodically. We design the seismic array by placing stations at the vertex or inside the tiles. To ensure a fair comparison, the tiles are scaled such that the array aperture and the number of stations are the same as the regular (Fig. 10b) and random (Fig. 15a) one. A total of five arrays are considered: The Hat and Specter array, where stations are placed at the vertex of the Hat and Specter tile; the central Hat array, where the stations are placed at the centre of the Hat tile; the Rectangular jittered array, where the stations are randomly jittered sampling of a rectangular array; as well as the Sparsely random array formed by the half down sampling of the Rectangular jittered array. Figs 19 and 20 show the station distribution and ARF of these five arrays. The same synthetic tests as Section 5.1 are used to assess the array performance in generating the quality of dispersion images. Fig. 21 shows the dispersion images generated using the original and new imaging conditions.

Figure 19.

The station (black triangles) distribution of the array designed by the tiles of ‘the Hat’ family. (a) The shape of the Hat tile. (b) The shape of the Turtle tile. (c)The shape of the Specter tile. (d) The station distribution and ARF of the Hat array. (e) The station distribution and ARF of the Specter array. (f) The station distribution and ARF of the central Hat array.

Open in new tab Download slide

Figure 20.

The station (black triangles) distribution and ARF of the random array. (a) The station distribution and ARF of the Rectangular jittered array. (b) The Sparsely random array, where the stations are the half down sampling of the Rectangular jittered array. The blue dots in left panels are the centre of the grids.

Open in new tab Download slide

Figure 21.

The dispersion images generated by different BF imaging conditions based on the synthetic data for the array shown in Figs 19 and 20. The left two columns show the results obtained by the original imaging conditions, while the right two columns show the results obtained by the new imaging conditions. For each imaging condition, the dispersion images are plotted in the f–k and f–v domain, respectively.

Open in new tab Download slide

Fig. 19(d) shows the largest side lobe amplitude is comparable to that of the main one for the Hat array (see also Fig. C1 in Appendix C for quantitative comparison). Fig. 21(a) shows the dispersion image for this array is characterized similarly to those for the regular array shown in Fig. 13 (only the positive values are colour-coded in Fig. 21). The Hat array does not perform as well as the random array. The artifacts associated with |${k_n}( \omega ) + m{k_{\textrm{aliasing}}}$| remain. This is not a surprise. As shown in Fig. 19(a), the vertices of the Hat tile are located at the centre of the tiling hexagon or the middle point of its side. This distribution results in a sampling pattern close to the regular array (see Fig. C2 in the Appendix C). Only the wavefield at few specific interstation distances are sampled. A similar feature can be observed if we put the stations on the vertex of the Turtle tile, which is also located at the specific position of the tiling hexagon (Fig 19b).

For the Specter array, the largest amplitude of the side lobe is less than half that of the main one (see also Fig. C1 in Appendix C). Fig. 21(b) shows that the dispersion image for this array is better than the Hat array. However, significant artifacts still remain due to the fact that some vertices of the Specter tile are located at the centre or the side midpoint of the tiling hexagon, resulting the sampling on the spatial wavefields exhibit weak periodicity (see Fig. C2 in Appendix C).

Fig. 19(f) shows that the ARF of the central Hat array outperforms the Hat and Specter arrays, and correspondingly a better dispersion image is expected, as shown in Fig. 21(c). However, the artifacts still remain.

In terms of the removal of artifacts associated with the regular sampling pattern, if the wavefield sampling becomes more non-uniform, the artifacts would be more diffuse (e.g. Zwartjes & Sacchi 2007) and incoherent. Although aperiodic tiling provides a deterministic idea for seismic array design, aperiodic tiling does not imply the sampling pattern of the spatial wavefield is non-periodic. For the spatial wavefield, it is sufficient to sample once at the same location, that is, the station pair with the same spacing appear once. This is obviously not the case of the aperiodic tiling. Therefore, the array performance designed based on the aperiodic monotile usually do not perform as well as the random array shown in Fig. 15(a).

To design the array randomly may provide an alternative, as illustrated in Section 6.1. In practice, however, we usually need to provide defined ‘random’ array layout for field worker. It might be a solution to generate a random array by jittering a regular array. Jittered sampling is introduced by Hennenfent & Herrmann (2008) for the 1-D case in seismic survey and is extended to the 2-D case by Shahidi et al. (2013). Jittered sampling limits the perturbations to a certain region, while using the randomness of the sampling to turn the aliasing into incoherent noise.

In Fig. 20(a), we use jittered sampling of the rectangular tiling to generate a random array. The blue circle indicates the centre of the grid. The random perturbation around the centre is restricted in the interior of the grid to ensure a relatively uniform distribution. Fig. 21(e) shows the results of the synthetic test, which, as also indicated by the ARF, is comparable to the result of the random array shown in Fig. 17.

Random down sampling is also used in exploration geophysics to reconstruct the wavefield using less information (e.g. Herrmann 2010). In Fig. 20(b), we down sample Fig. 20(a) by randomly retaining half of the stations. Fig. 21(e) shows the corresponding results of the synthetic test. With half of stations, we obtain almost the same results as in Fig. 21(d).

According to the theory of compressed sensing (e.g. Donoho 2006; Baraniuk 2007), the complete seismic wavefield can be reconstructed based on random sparse sampling, which is not restricted by the Nyquist sampling criterion. The prerequisite is that the wavefield is sparse in a proper domain. As a function of frequency, the phase velocities of seismic surface waves are eigenvalues of the surface wave dispersion equation, and they exist only at specific wavenumbers and frequencies. The surface wave mode branches satisfy sparsity in the f–k domain. In the framework of compressed sensing theory, we may only need fewer random samples to reconstruct the surface wave wavefield without aliasing. The application of compressed sensing to the current topic is beyond the scope of this paper. However, the above investigation shows that random arrays outperform the regular arrays and the one designed based on the Hat family in aperiodic tiling. Jittered sampling of a regular array is a feasible design option for random array.

7 REMOVAL OF RADIAL ARTIFACTS

Although the radial artifacts associated with |$m{k_{\textrm{aliasing}}}$| is also related to the array configuration, its appearance depends on the wavefield characteristics. It would be observed if the amplitude related to the constant wavenumber is comparable to the surface wave of interest. In this section, we show they can be reduced by pre-processing the NCFs.

As shown in Fig. 22(a), the array is a subset of ChinArray (Phase II) containing 104 stations (blue triangles) with an average interval of 30 km. Wang et al. (2020) have extracted the fundamental mode Rayleigh wave between 7 and 35 s using the conventional CBF. Based on the NCFs, which are the stacking results of three months continuous records (from 2014 January 1 to 2014 March 31), we investigate the extraction of higher modes by MCBF. Fig. 22(b) shows the NCFs with bandpass filter of 0.01–0.5 Hz. The waveforms with much larger amplitude are observed around zero-lag time. Figs 22(c) and (d) show the dispersion image generated by eqs (27) and (26), respectively. The white dashed lines denote the radial artifacts predicted using ARF. Although the crossed artifacts have been eliminated in Fig. 22(d) with the new imaging conditions, the improvement of the images is not apparent due to the fact that the radial artifacts are much stronger than the crossed one.

Figure 22.

An illustration of the radial artifacts of the dispersion image in f–v domain. (a) The subset (blue triangles) of ChinArray (Phase II). (b) The interstation NCFs. (c) The dispersion image generated by original imaging condition. (d) The dispersion image generated by new imaging condition. The white dashed lines denote the radial artifacts predicted by ARF. New imaging condition is effective for the removal of crossed artifacts but not for radial artifacts. The black dashed lines represent the dispersion curves predicted by the averaged model given by Wang et al. (2020).

Open in new tab Download slide

To remove the radial artifacts, we eliminate the waveforms near the zero-lag time and keep those with velocity lower than 5.5 km s⁻¹, as shown in Fig. 23(a). The corresponding dispersion image is given in Fig. 23(b), where the radial artifacts are greatly reduced. We then use a cosine-taper window to select the waveform of surface waves and suppress the disturbance at higher frequencies. Fig. 23(c) shows the windowed NCFs using the velocity range from 2 to 5.5 km s⁻¹. The corresponding dispersion image is given in Fig. 23(d). Compared with Fig. 22(b), the disturbances at high frequency are significantly reduced. The windowing operation reduces radial artifacts and improves the quality of dispersion image, especially when a significant waveform around zero lag time are observed in NCFs.

An illustration of the removal of radial artifacts. (a) The interstation NCFs after erasing the waveform with group velocity greater than 5.5 km s−1. (b) The dispersion image generated by new imaging condition based on the NCFs in panel (a). (c) The time windowed NCFs after keeping only the waveforms with group velocity from 2 to 5.5 km s−1. (d) The same as (b) but based on the NCFs in panel (c). The white dashed lines are the predicted artifacts.

Figure 23.

An illustration of the removal of radial artifacts. (a) The interstation NCFs after erasing the waveform with group velocity greater than 5.5 km s⁻¹. (b) The dispersion image generated by new imaging condition based on the NCFs in panel (a). (c) The time windowed NCFs after keeping only the waveforms with group velocity from 2 to 5.5 km s⁻¹. (d) The same as (b) but based on the NCFs in panel (c). The white dashed lines are the predicted artifacts.

Open in new tab Download slide

The black dashed lines in Figs 22(d) and 23(d) represent the dispersion curves predicted by the averaged model given by Wang et al. (2020). For the fundamental mode Rayleigh wave, the predicted dispersion curve agrees well with the observed one in the frequency range of 0.03–0.15 Hz. This is not surprised since the model in Wang et al. (2020) is derived from the inversion of the fundamental mode Rayleigh dispersion curve at this frequency range. For the fundamental mode Rayleigh wave at frequencies higher than 0.15 Hz and the higher modes, the predicted dispersion curve deviated from the observed one. This also implies the importance of high modes for constraint of the deep structure.

8 CONCLUSIONS AND DISCUSSIONS

In this paper, we study how to reduce the artifacts of the multimode dispersion image generated by beamforming the ambient seismic noise. Our investigation into the resolution and aliasing effects reveals some insights. The main contributions are:

We identify three types of artifacts in the multimode dispersion image obtained by beamforming, namely, the artifacts associated with periodic extension of the positive wavenumber, the crossed artifacts associated with the negative wavenumber and those radial artifacts associated with constant wavenumber. Formulas for predicting them are provided.
New BF imaging conditions are proposed to eliminate crossed artifacts. Compared with the existing methods, the new imaging conditions directly utilize the real and imaginary parts of the NCFs without constructing the complete GF via the Hilbert transform.
We investigate the array performance of random array and the array designed by the tiles of the Hat family, a newly discovery in the field of monotile aperiodic tiling, showing that the random array would help to reduce the artifacts associated with the periodic extension of the negative wavenumber, as well as the positive wavenumber. Random arrays outperform the regular arrays and the one designed based on the Hat family in aperiodic tiling. Jittered sampling of a regular array is a feasible design option for random array.
We show that adding a window to the NCFs eliminates the radial artifacts associated with the constant wavenumber.

Thanks to scientific progress and economic development, dense array observations and array-based high-resolution techniques have been popularized in the field of seismic imaging. Dense array observation duration is usually insufficient to record enough seismic events, the noise-based method is hence becoming general for array-based technique, especially since it has been found the high-mode surface wave can be extracted by array-based method from noise recordings. Regular array suffers from aliasing for array-based surface wave imaging. Aperiodic tiling provides a scheme for array design, but it does not perform as well as random arrays as far as artifact removal is concerned.

Random array design seems to be a paradox. On the one hand, due to the complex environmental and geological conditions, it is usually impossible to lay out the stations according to the pre-designed regular array in the field work. The advantages of random array free us from the worry of not being able to implement a designed regular pattern. On the other hand, from the viewpoint for providing practical guidance on the deployment of the array, we need at least a defined ‘random’ array pattern for field workers. Two suggested principles for the design of random arrays are: (1) The sampling density should be as uniform as possible. (2) Sampling should not be repetitive or regular whenever possible. A practical design scheme would be a randomly jittered based on a uniform regular array, as suggested in this paper.

Furthermore, the discussion of array performance and random array design in this paper is specific to individual array. For array-based imaging, we first divide a large array into many subarrays and then combine the results from each subarray to image the structure under the entire array. Similar to the effect imposed by the imprints of the ray coverage on the result obtained by ray-based method, the array-based imaging results may also be smeared by the subarray division and geometry. When designing a random array, it would be of great value to consider simultaneously the performance of individual subarray, the station distribution of the entire array, as well as the subarray division in the imaging scheme.

ACKNOWLEDGEMENTS

We thank the editors, Michal Malinowski and Christoph Sens-Schönfelder, and the reviewers, Feng Cheng, Roméo Courbis and an anonymous one, for their constructive comments and suggestions, which helped the authors to improve the manuscript greatly. This work was supported by the National Natural Science Foundation of China (42474081), the Fundamental Research Funds of the Institute of Geophysics, China Earthquake Administration (DQJB23K34), and the National Key R & D Program of China (2017YFC1500200).

DATA AVAILABILITY

The NCFs of the array used as the examples in this paper, as well as the Python package for generating the dispersion image using new BF imaging conditions based on these NCFs, are openly available from 10.6084/m9.figshare.25678170.v1 (Qin 2024).

REFERENCES

Abramowitz

Stegun

I.A.

1964

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, NBS Applied Mathematics Series 55

National Bureau of Standards

Washington, DC

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Aki

1957

Space and time spectra of stationary stochastic waves, with special reference to microtremors

Bull. Earthq. Res. Inst.

415

–

456

Google Scholar

OpenURL Placeholder Text

WorldCat

Asten

M.W.

Henstridge

J.D.

1984

Array estimators and the use of microseisms for reconnaissance of sedimentary basins

Geophysics

(

1828

–

1837

Baraniuk

R.G.

2007

Compressive sensing

IEEE Signal Process. Mag.

(

118

–

121

10.1109/MSP.2007.4286571

Bensen

G.D.

et al. ,

2007

Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements

Geophys. J. Int.

169

(

1239

–

1260

10.1111/j.1365-246X.2007.03374.x

Boschi

Weemstra

Verbeke

Ekström

Zunino

Giardini

2013

On measuring surface wave phase velocity from station–station cross-correlation of ambient signal

Geophys. J. Int.

192

(

346

–

358

Campillo

Paul

2003

Long-range correlations in the diffuse seismic coda

Science

299

(

5606

547

–

549

10.1126/science.1078551

Capon

1969

High-resolution frequency-wavenumber spectrum analysis

Proc. IEEE

(

1408

–

1418

10.1109/PROC.1969.7278

Chen

Liu

Cheng

2023

Multi-mode imaging of ambient background noise for karst detection in the limestone area based on frequency-bessel transform

Appl. Sci.

(

5135

Cheng

Xia

2023

Artifacts in high-frequency passive surface wave dispersion imaging: toward the linear receiver array

Surv. Geophys.

(

1009

–

1039

10.1007/s10712-023-09772-1

Cheng

Xia

Zhang

Zhou

Ajo-Franklin

J.B.

2021

Phase-weighted slant stacking for surface wave dispersion measurement

Geophys. J. Int.

226

(

256

–

269

Donoho

D.L.

2006

Compressed sensing

IEEE Trans. Inf. Theory

(

1289

–

1306

10.1109/TIT.2006.871582

Forbriger

2003

Inversion of shallow-seismic wavefields: I. Wavefield transformation

Geophys. J. Int.

153

(

719

–

734

10.1046/j.1365-246X.2003.01929.x

Foti

et al.

2018

Guidelines for the good practice of surface wave analysis: a product of the InterPACIFIC project

Bull. Earthq. Eng.

(

2367

–

2420

10.1007/s10518-017-0206-7

Pan

Dong

Chen

2022

Improved high-resolution 3D vs model of Long Beach, CA: inversion of multimodal dispersion curves from ambient noise of a dense array

Geophys. Res. Lett.

(

e2021GL097619

Gal

Reading

A.M.

Ellingsen

S.P.

2019

Short timescale analysis of microseisms and application to array calibration

J. Geophys. Res. Solid Earth

124

(

2684

–

2701

Gerstoft

Tanimoto

2007

A year of microseisms in southern California

Geophys. Res. Lett.

(

2007GL031091

10/d2sts6

Gradshteyn

I.S.

Ryzhik

I.M.

2007

Table of Integrals, Series, and Products

Academic Press, Elsevier Inc

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Haney

M.M.

Mikesell

T.D.

Van Kasper

2012

Extension of the spatial autocorrelation (SPAC) method to mixed-component correlations of surface waves

Geophys. J. Int.

191

189

–

206

10.1111/j.1365-246X.2012.05597.x

Harmon

Gerstoft

Rychert

C.A.

Abers

G.A.

Cruz

M.S.

Fischer

K.M.

2008

Phase velocities from seismic noise using beamforming and cross-correlation in Costa Rica and Nicaragua

Geophys. Res. Lett.

(

2008GL035387

10/cjgfrp

Hennenfent

Herrmann

F.J.

2008

Simply denoise: wavefield reconstruction via jittered undersampling

Geophysics

(

V19

–

V28

Herrmann

F.J.

2010

Randomized sampling and sparsity: getting more information from fewer samples

Geophysics

(

173

–

187

Horike

1985

Inversion of phase velocity of long-period microtremors to the S-wave-velocity structure down to the basement in urbanized areas

J. Phys. Earth

(

–

10.4294/jpe1952.33.59

Luo

Yao

2020

The frequency-bessel spectrograms of multicomponent cross-correlation functions from seismic ambient noise

J. Geophys. Res. Solid Earth

125

(

e2020JB019630

Lacoss

R.T.

Kelly

E.J.

Toksöz

M.N.

1969

Estimation of seismic noise structure using arrays

Geophysics

(

–

Shi

Chen

2021

Constraints on crustal P wave structure with leaking mode dispersion curves

Geophys. Res. Lett.

(

e2020GL091782

Liu

Wang

Chang

Zhu

2023

Rayleigh wave phase velocity maps at regional scale inferring from SPAC of ambient noise at a dense array: a case study in northeastern Tibetan plateau

Pure Appl. Geophys.

180

(

1973

–

1988

10.1007/s00024-023-03266-2

Lobkis

O.I.

Weaver

R.L.

2001

On the emergence of the Green's function in the correlations of a diffuse field

J. Acoust. Soc. Am.

110

(

3011

–

3017

Löer

Riahi

Saenger

E.H.

2018

Three-component ambient noise beamforming in the Parkfield area

Geophys. J. Int.

213

(

1478

–

1491

Löer

Toledo

Norini

Zhang

Curtis

Saenger

E.H.

2020

Imaging the deep structures of Los Humeros geothermal field, Mexico, using three-component seismic noise beamforming

Seismol. Res. Lett.

(

3269

–

3277

2021

Revisiting the cross-correlation and SPatial AutoCorrelation (SPAC) of the seismic ambient noise based on the plane wave model

Rev. Geophys. Planet. Phys.

(

123

–

163

10.19975/j.dqyxx.2020-016

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Wang

Ding

2018

Effect of uneven noise source and/or station distribution on estimating the azimuth anisotropy of surface waves

Earthq. Sci.

(

175

–

186

10.29382/eqs-2018-0175-1

Luo

Zhou

Yao

2022

Improvement of frequency–Bessel phase-velocity spectra of multicomponent cross-correlation functions from seismic ambient noise

Bull. seism. Soc. Am.

112

(

2257

–

2279

Mordret

Grushin

A.G.

2024

Beating the aliasing limit with aperiodic monotile arrays

. https://arxiv.org/abs/2408.16476v1,

(accessed December 01, 2024).

OpenURL Placeholder Text

WorldCat

Nakahara

2006

A systematic study of theoretical relations between spatial correlation and Green's function in one-, two- and three-dimensional random scalar wavefields

Geophys. J. Int.

167

(

1097

–

1105

10.1111/j.1365-246X.2006.03170.x

Nakata

Gualtieri

Fichtner

2019

Seismic ambient noise

, Cambridge:

Cambridge University Press

Näsholm

S.P.

Iranpour

Wuestefeld

Dando

B.D.E.

Baird

A.F.

Oye

2022

Array signal processing on distributed acoustic sensing data: directivity effects in slowness space

J. Geophys. Res. Solid Earth

127

(

e2021JB023587

Nimiya

Ikeda

Tsuji

2023

Multimodal rayleigh and Love wave joint inversion for S-wave velocity structures in Kanto Basin, Japan

J. Geophys. Res. Solid Earth

128

e2022JB025017

O'Neill

M.E.

2014

PCG: a family of simple fast space-efficient statistically good algorithms for random number generation (No. HMC-CS-2014-0905)

Harvey Mudd College

. https://www.cs.hmc.edu/tr/hmc-cs-2014-0905.pdf,

(accessed February 6 2025)

Picozzi

Parolai

Bindi

2010

Deblurring of frequency-wavenumber images from small-scale seismic arrays

Geophys. J. Int.

181

(

357

–

368

10.1111/j.1365-246X.2009.04471.x

Qin

2024

. NCFs, python package and sample script of new imaging condition for beamforming,

10.6084/m9.figshare.25678170.v1

[Dataset]

(accessed February 6 2025)

Qin

2024

Improved beamforming schemes for estimation of multimode surface wave dispersion curves from seismic noise with reducing effect of the irregular array geometry and/or anisotropic source distribution

Geophys. J. Int.

237

(

250

–

270

Qin

Ding

Feng

Zhang

2022

High-resolution 3D shallow S wave velocity structure of Tongzhou, subcenter of Beijing, inferred from multi-mode rayleigh waves by beamforming seismic noise at a dense array

J. Geophys. Res. Solid Earth

127

(

e2021JB023689

Riahi

Saenger

E.H.

2014

Rayleigh and Love wave anisotropy in Southern California using seismic noise

Geophys. Res. Lett.

363

–

369

Rost

Thomas

2002

Array seismology: methods and applications

Reviews of Geophysics

(

1008

10.1029/2000RG000100

Roux

Ben-Zion

2017

Rayleigh phase velocities in Southern California from beamforming short-duration ambient noise

Geophys. J. Int.

211

(

450

–

454

Ruigrok

Gibbons

Wapenaar

2017

Cross-correlation beamforming

J. Seismol.

(

495

–

508

10.1007/s10950-016-9612-6

Seydoux

de Rosny

Shapiro

N.M.

2017

Pre-processing ambient noise cross-correlations with equalizing the covariance matrix eigenspectrum

Geophys. J. Int.

210

(

1432

–

1449

Shahidi

Tang

Herrmann

F.J.

2013

Application of randomized sampling schemes to curvelet-based sparsity-promoting seismic data recovery

Geophys. Prospect.

973

–

997

10.1111/1365-2478.12050

Smith

Myers

J.S.

Kaplan

C.S.

Goodman-Strauss

2024

A chiral aperiodic monotile

10.5070/C64264241

(accessed February 6 2005

OpenURL Placeholder Text

WorldCat

Crossref

Snieder

Wakin

2022

When randomness helps in undersampling

SIAM Rev.

(

1062

–

1080

Takagi

Nishida

2022

Multimode dispersion measurement of surface waves extracted by multicomponent ambient noise cross-correlation functions

Geophys. J. Int.

231

(

1196

–

1220

Tokimatsu

1997

Geotechnical site characterization using surface waves

Proc. 1st Intl. Conf. Earthquake Geotechnical Engineering

1333

–

1368

Ishihara

A.A. Balkema, Rotterdam

Netherlands

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Tsai

V.C.

2009

On establishing the accuracy of noise tomography travel-time measurements in a realistic medium

Geophys. J. Int.

178

(

1555

–

1564

10.1111/j.1365-246X.2009.04239.x

Tsai

V.C.

Moschetti

M.P.

2010

An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results

Geophys. J. Int.

182

(

454

–

460

10.1111/j.1365-246X.2010.04633.x

Tsarsitalidou

Hillers

Giammarinaro

Boué

Stehly

Campillo

2024

Long period rayleigh wave focal spot imaging applied to USArray data

J. Geophys. Res. Solid Earth

129

(

e2023JB027417

Viens

Perton

Spica

Z.J.

Nishida

Yamada

Shinohara

2023

Understanding surface wave modal content for high-resolution imaging of submarine sediments with distributed acoustic sensing

Geophys. J. Int.

232

(

1668

–

1683

Wang

Chen

2019

Frequency-bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data

J. Geophys. Res. Solid Earth

124

(

3708

–

3723

Wang

Maupin

Ding

Zheng

Zhong

2020

Surface wave tomography of northeastern Tibetan plateau using beamforming of seismic noise at a dense array

J. Geophys. Res. Solid Earth

125

(

e2019JB018416

Wapenaar

Fokkema

2006

Green's function representations for seismic interferometry

Geophysics

(

SI33

–

SI46

Wathelet

Jongmans

Ohrnberger

Bonnefoy-Claudet

2008

Array performances for ambient vibrations on a shallow structure and consequences over vs inversion

J. Seismol.

(

–

10.1007/s10950-007-9067-x

Woods

J.W.

Lintz

P.R.

1973

Plane waves at small arrays

Geophysics

(

1023

–

1041

Wang

2023a

Ambient noise surface wave tomography of quaternary structures derived from a high-density array in the central Hebei Depression, North China

Geosci. J.

(

177

–

190

10.1007/s12303-022-0033-y

Wang

2023b

Multi-mode surface wave tomography of a water-rich layer of the Jizhong depression using beamforming at a dense array

Remote. Sens.

(

Article 1

Xia

Dai

Liu

Ning

2021

Modified frequency–Bessel transform method for dispersion imaging of Rayleigh waves from ambient seismic noise

Geophys. J. Int.

225

(

1271

–

1280

Yamaya

Mochizuki

Akuhara

Nishida

2021

Sedimentary structure derived from multi-mode ambient noise tomography with dense OBS network at the Japan trench

J. Geophys. Res. Solid Earth

126

(

e2021JB021789

Yao

Cao

Huang

2023

Automatic extraction of surface wave dispersion curves using unsupervised learning and high-resolution tau-p transform

Earth Space Sci.

(

e2023EA003198

Yokoi

Margaryan

2008

Consistency of the spatial autocorrelation method with seismic interferometry and its consequence

Geophys. Prospect.

435

–

451

10.1111/j.1365-2478.2008.00709.x

Zhou

Chen

2022

Removal of crossed artifacts from multimodal dispersion curves with modified frequency–Bessel method

Bull. seism. Soc. Am.

112

(

143

–

152

Zwartjes

P.M.

Sacchi

M.D.

2007

Fourier reconstruction of nonuniformly sampled, aliased seismic data

Geophysics

(

V21

–

V32

APPENDIX A: CONVENTION FOR THE FOURIER TRANSFORM, HILBERT TRANSFORM AND CROSS-CORRELATION

The following convention for the Fourier Transform from the time (t) domain to the frequency (⁠|$\omega $|⁠) domain is adopted in this paper

$$\begin{eqnarray} \begin{array}{@{}l@{}} f(t) = \int\limits_{ - \infty }^\infty {F(\omega ){e^{\iota \omega t}}\mathrm{ d}\omega } \\ F(\omega ) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^\infty {f(t){e^{ - \iota \omega t}}\mathrm{ d}t} \end{array} . \end{eqnarray}$$

(A1)

Accordingly, for the Fourier transform from the space (r) domain to the wavenumber domain (k) we use the convention

$$\begin{eqnarray} \begin{array}{@{}l@{}} f(r) = \int\limits_{ - \infty }^\infty {F(k){e^{ - \iota kr}}\mathrm{ d}k} \\ F(k) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^\infty {f(r){e^{\iota kr}}\mathrm{ d}r} \end{array} \end{eqnarray}$$

(A2)

|$f(t)$| and |$F(\omega )$|⁠, as well as |$f(r)$| and |$F(k)$| are two Fourier transform pairs. These conventions imply the plane wave expression |${e^{ - ikx}}{e^{i\omega t}}$| describes the wave propagating in the positive x, while |${e^{ + ikx}}{e^{i\omega t}}$| describes the wave propagating in the negative x. For the cylindrical wave in the Cartesian coordinate system with z downwards, the zero-order Hankel function |$H_{\rm{0}}^{{\rm{(1)}}}(kr)$| of the first kind represents the converging wave and the zero-order Hankel function |$H_{\rm{0}}^{{\rm{(2)}}}(kr)$| of the second kind represents the wave propagating outwards.

The Hilbert transform of the function |$s(x)$| is defined as

$$\begin{eqnarray} s(t) = \mathcal{H}\ [s(x)] = \frac{1}{\pi }{\rm{P}}{\rm{.V}}.\int\limits_{ - \infty }^{ + \infty } {\frac{{s(x)}}{{t - x}}\mathrm{ d}\tau } , \end{eqnarray}$$

(A3)

where P.V. represents the principal value of the Cauchy integral. Under this definition, the analytical signal |$S(t)$|of the real-value signal |$s(t)$| can be expressed as

$$\begin{eqnarray} S(t) = s(t) + \iota \mathcal{H} [s(t)] \end{eqnarray}$$

(A4)

The symbol |$\mathcal{H}$|represents the Hilbert transform.

The cross-correlation of two (complex) signals is defined as:

$$\begin{eqnarray} {C_{fg}}(\tau ) = \int\limits_{ - \infty }^{ + \infty } {{f^{\rm{*}}}(t)g(t + \tau )\mathrm{ d}t} , \end{eqnarray}$$

(A5)

where superscript * represents complex conjugate. If the Fourier transform of |$g(t)$| is |$G(\omega )$|⁠. According to this definition and the above convention for Fourier transform, we have

$$\begin{eqnarray} {C_{fg}}(\tau ) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^{ + \infty } {{F^*}(\omega )} G(\omega ){e^{i\omega \tau }}\mathrm{ d}\omega \end{eqnarray}$$

(A6)

i.e. |${C_{fg}}(\tau )$| and |${F^*}(\omega )G(\omega )$| are Fourier transform pair. Cross-correlation in the time domain corresponds to the product in the frequency domain by taking complex conjugate of one of them.

APPENDIX B: THEORETICAL REPRESENTATION OF MCBF

For an ideal cylindrical wavefield, substituting |$C(r,\omega ) = {J_0}({k_n}r)$| into eq. (4). The azimuth-averaged MCBF can be written as

$$\begin{eqnarray} \overline {\textrm{MCBF}} \left( {k,\omega } \right) &=& \frac{1}{{{N^2}}}\sum\limits_i^N {\sum\limits_j^N {\sqrt {\pi k{r_{ij}}} {J_0}\left( {{k_n}{r_{ij}}} \right){e^{\iota \left( {k{r_{ij}} - \frac{\pi }{4}} \right)}}} } \\ &=& \frac{{\sqrt 2 }}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( {\sum\limits_i^N {\sum\limits_j^N {\cos \left( {{k_n}{r_{ij}} - \frac{\pi }{4}} \right){e^{\iota \left( {k{r_{ij}} - \frac{\pi }{4}} \right)}}} } } \right) . \end{eqnarray}$$

(B1)

In eq. (B1), the far field approximation of |${J_0}({k_n}r)$| is applied. The double summation over the number of stations means the summation over the interstation distance |${r_{ij}}$|⁠. Since the geometric spread is corrected by multiplying |${r_{ij}}$| in eq. (B1), the autocorrelation for |$i = j$| is excluded. The cross-correlations for |$i \ne j$| are counted twice for the same interstation distance |${r_{ij}} = {r_{ji}}$|⁠. Eq. (B1) can then be written as

$$\begin{eqnarray} \overline {\textrm{MCBF}} \left( {k,\omega } \right) &=& \frac{{2\sqrt 2 }}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\sum\limits_j^M {\cos \left({k_n}{r_j} - \frac{\pi }{4}\right){e^{\iota \left(k{r_j} - \frac{\pi }{4}\right)}}} \\ & =& \frac{{2\sqrt 2 }}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\sum\limits_j^M {\cos \left({k_n}{r_j} - \frac{\pi }{4}\right)\left( {\cos \left(k{r_j} - \frac{\pi }{4}\right) + \iota \sin \left(k{r_j} - \frac{\pi }{4}\right)} \right)} \\ & =& \mathrm{ CC} + \iota \mathrm{ CS} , \end{eqnarray}$$

(B2)

where |$M = {{N(N - 1)} {/ {\vphantom {{N(N - 1)} 2}} } 2}$|⁠. Let

$$\begin{eqnarray} \mathrm{ CC}^{\prime} = \sum\limits_j^M {\cos \left({k_n}{r_j} - \frac{\pi }{4}\right)\cos \left(k{r_j} - \frac{\pi }{4}\right)} = \frac{1}{2}\sum\limits_j^M {\left[ {\cos ({k_n}{r_j} - k{r_j}) + \sin ({k_n}{r_j} + k{r_j})} \right]} . \end{eqnarray}$$

(B3)

It is assumed the cylindrical wavefield is sampled spatially in radial distance with equal interval |$\Delta r$|⁠. We have |${r_j} = j\Delta r$|⁠. Eq. (B3) can be rewritten as:

$$\begin{eqnarray} \mathrm{ CC}^{\prime} = \frac{1}{2}\sum\limits_j^M {\left[ {\cos ({k_n} - k)j\Delta r + \sin ({k_n} + k)j\Delta r} \right]} . \end{eqnarray}$$

(B4)

Applying the eqs (1.342.1) and (1.342.2) in Gradshteyn & Ryzhik (2007)

$$\begin{eqnarray} \sum\limits_{k = 0}^n {\cos kx} &=& \cos {\textstyle{{n + 1} \over 2}}x\frac{{\sin {\textstyle{{nx} \over 2}}}}{{\sin {\textstyle{x \over 2}}}} + 1 \\ \sum\limits_{k = 1}^n {\sin kx} &=& \sin {\textstyle{{n + 1} \over 2}}x\frac{{\sin {\textstyle{{nx} \over 2}}}}{{\sin {\textstyle{x \over 2}}}} . \end{eqnarray}$$

(B5)

Eq. B4 can be recast into

$$\begin{eqnarray} \mathrm{ CC}^{\prime} &=& \frac{1}{2}\left[ {\cos {\textstyle{{(M + 1)({k_n} - k)\Delta r} \over 2}}\frac{{\sin {\textstyle{{M({k_n} - k)\Delta r} \over 2}}}}{{\sin {\textstyle{{({k_n} - k)\Delta r} \over 2}}}} + \sin {\textstyle{{(M + 1)({k_n} + k)\Delta r} \over 2}}\frac{{\sin {\textstyle{{M({k_n} + k)\Delta r} \over 2}}}}{{\sin {\textstyle{{({k_n} + k)\Delta r} \over 2}}}}} \right] \\ & =& \frac{M}{2}\left[ {\cos {\textstyle{{(M + 1)({k_n} - k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} - k)\Delta r} \over 2}}} \right) + \sin {\textstyle{{(M + 1)({k_n} + k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} + k)\Delta r} \over 2}}} \right)} \right] , \end{eqnarray}$$

(B6)

where

$$\begin{eqnarray} {D_M}(x) = \frac{{\sin {\textstyle{{Mx} \over 2}}}}{{M\sin {\textstyle{x \over 2}}}} \end{eqnarray}$$

(B7)

$$\begin{eqnarray} \mathrm{ CC}^{\prime} \approx \frac{M}{2}\left[ {\cos {\textstyle{{R({k_n} - k)} \over 2}}{\rm{sinc}}{\textstyle{{R({k_n} - k)} \over 2}} + \sin {\textstyle{{R({k_n} + k)} \over 2}}{\rm{sinc}}\left( {{\textstyle{{R({k_n} + k)} \over 2}}} \right)} \right] , \end{eqnarray}$$

(B8)

where |${\rm{sinc(}}x) = {{\sin x} {/ {\vphantom {{\sin x} x}} } x}$| is the sinc function. Similarly, it can be deduced

$$\begin{eqnarray} \mathrm{ CS}^{\prime} &=& \sum\limits_j^M {\cos ({k_n}{r_j} - \frac{\pi }{4})\sin (k{r_j} - \frac{\pi }{4})} = \frac{1}{2}\sum\limits_j^M {\left[ {\cos ({k_n}{r_j} + k{r_j}) - \sin ({k_n}{r_j} - k{r_j})} \right]} \\ &=& \frac{M}{2}\left[ {\cos {\textstyle{{(M + 1)({k_n} + k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} + k)\Delta r} \over 2}}} \right) - \sin {\textstyle{{(M + 1)({k_n} - k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} - k)\Delta r} \over 2}}} \right)} \right] \\ & \approx &\frac{M}{2}\left[ {\cos {\textstyle{{R({k_n} + k)} \over 2}}{\rm{sinc}}{\textstyle{{R({k_n} + k)} \over 2}} - \sin {\textstyle{{R({k_n} - k)} \over 2}}{\rm{sinc}}\left( {{\textstyle{{R({k_n} - k)} \over 2}}} \right)} \right] , \end{eqnarray}$$

(B9)

$$\begin{eqnarray} \mathrm{ SS}^{\prime} &=& \sum\limits_j^M {\sin ({k_n}{r_j} - \frac{\pi }{4})\sin (k{r_j} - \frac{\pi }{4})} = \frac{1}{2}\sum\limits_j^M {\left[ {\cos ({k_n}{r_j} - k{r_j}) - \sin ({k_n}{r_j} + k{r_j})} \right]} \\ & =& \frac{M}{2}\left[ {\cos {\textstyle{{(M + 1)({k_n} - k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} - k)\Delta r} \over 2}}} \right) - \sin {\textstyle{{(M + 1)({k_n} + k)\Delta r} \over 2}}{D_M}\left( {{\textstyle{{({k_n} + k)\Delta r} \over 2}}} \right)} \right] \\ & \approx & \frac{M}{2}\left[ {\cos {\textstyle{{R({k_n} - k)} \over 2}}{\rm{sinc}}{\textstyle{{R({k_n} - k)} \over 2}} - \sin {\textstyle{{R({k_n} + k)} \over 2}}{\rm{sinc}}\left( {{\textstyle{{R({k_n} + k)} \over 2}}} \right)} \right] . \end{eqnarray}$$

(B10)

The modulus of |$\mathrm{ CC} + \iota \mathrm{ CS}$| can be approximated as

$$\begin{eqnarray} \left| {\mathrm{ CC} + \iota \mathrm{ CS}} \right| = \frac{{\sqrt 2 }}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( {{\rm{sinc}}\frac{{R({k_n} + k)}}{2} + {\rm{sinc}}\frac{{R({k_n} - k)}}{2}} \right) . \end{eqnarray}$$

(B11)

The theoretical representations of the azimuth-averaged MCBF with different imaging conditions can then be written as

$$\begin{eqnarray} \mathrm{ CC} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( {\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} + \sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}} \right) \\ \mathrm{ CS} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( {\sin \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} - \cos \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}} \right) \\ \mathrm{ SC} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( { - \sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2} - \cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2}} \right) \\ \mathrm{ SS} &=& \frac{{\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\left( {\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} - \sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2}} \right) \\ \mathrm{ CC} + \mathrm{ SS} &=& \frac{{2\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\cos \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} \\ \mathrm{ CC} - \mathrm{ SS} &=& \frac{{2\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\sin \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2} \\ \mathrm{ CS} - \mathrm{ SC} &=& \frac{{2\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\sin \frac{{R({k_n} - k)}}{2}{\rm{sinc}}\frac{{R({k_n} - k)}}{2} \\ - \mathrm{ CS} - \mathrm{ SC} &=& \frac{{2\sqrt 2 M}}{{{N^2}}}\frac{{\sqrt k }}{{\sqrt {{k_n}} }}\cos \frac{{R({k_n} + k)}}{2}{\rm{sinc}}\frac{{R({k_n} + k)}}{2} . \end{eqnarray}$$

(B12)

APPENDIX C: ARF AND WAVEFIELD SAMPLING

Fig. C1 shows the cross-sections of the ARF of the arrays shown in Figs 19 and 20, where the ARF are plotted as a function of radial distance every two degrees. Compared to those 2-D displayed in Figs 19 and 20, the amplitude difference of the main and side lobes can be observed quantitatively along the cross-sections at different azimuth. This difference provides a criterion for evaluating array performance.

Figure C1.

The cross-sections of the ARF of the arrays shown in Figs 19 and 20. The grey lines denote the cross-sections of ARF, which are plotted as a function of radial distance every two degrees. The red dashed line denote the ARF after averaging over the azimuth.

Open in new tab Download slide

Fig. C2 shows the distribution of the sampling points for the array shown in Figs 19 and 20. To plot them, we first express the station-pairs with different orientations and interstation distances as the vectors in cartesian coordinate. A black dot is then plotted at the end of the vector to indicate the sampling point of that station-pair on the wavefield.

Figure C2.

The distribution of the sampling points for the array shown in Figs 19 and 20. For (a) and (b), only the points at the specific locations are sampled regularly. For (c) and (d), the sampling pattern appears with weak periodicity. For (e) and (f), the sampling points are randomly distributed with sufficient homogeneity.

Open in new tab Download slide

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Download all slides

Month:	Total Views:
February 2025	127
March 2025	205
April 2025	87

Article Contents

Artifacts removal from multimode dispersion image: new imaging conditions of beamforming and random array

SUMMARY

1 INTRODUCTION

2 BEAMFORMING

2.1 Beamforming schemes for extracting multimode dispersion curves

2.2 Dispersion image generated by different beamforming schemes

3 ARF, RESOLUTION AND ALIASING

3.1 ARF and resolution

3.2 Subpeaks and aliasing

3.3 Artifacts in the dispersion image

3.4 Removal of crossed artifacts based on Kramers–Kronig relations

4 NEW IMAGING CONDITIONS OF BEAMFORMING

4.1 Causal and acausal parts of GF

4.2 New imaging conditions of MCBF for removal of crossed artifacts

4.3 New imaging conditions of WCBF

5 EXAMPLES AND DISCUSSION

5.1 Synthetic data

5.2 Field data

6 RANDOM SAMPLING AND ARRAY DESIGN

6.1 Random sampling helps reduce aliasing: synthetic and field data

6.2 Random array design and performance comparison

7 REMOVAL OF RADIAL ARTIFACTS

8 CONCLUSIONS AND DISCUSSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: CONVENTION FOR THE FOURIER TRANSFORM, HILBERT TRANSFORM AND CROSS-CORRELATION

APPENDIX B: THEORETICAL REPRESENTATION OF MCBF

APPENDIX C: ARF AND WAVEFIELD SAMPLING

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Artifacts removal from multimode dispersion image: new imaging conditions of beamforming and random array

SUMMARY

1 INTRODUCTION

2 BEAMFORMING

2.1 Beamforming schemes for extracting multimode dispersion curves

2.2 Dispersion image generated by different beamforming schemes

3 ARF, RESOLUTION AND ALIASING

3.1 ARF and resolution

3.2 Subpeaks and aliasing

3.3 Artifacts in the dispersion image

3.4 Removal of crossed artifacts based on Kramers–Kronig relations

4 NEW IMAGING CONDITIONS OF BEAMFORMING

4.1 Causal and acausal parts of GF

4.2 New imaging conditions of MCBF for removal of crossed artifacts

4.3 New imaging conditions of WCBF

5 EXAMPLES AND DISCUSSION

5.1 Synthetic data

5.2 Field data

6 RANDOM SAMPLING AND ARRAY DESIGN

6.1 Random sampling helps reduce aliasing: synthetic and field data

6.2 Random array design and performance comparison

7 REMOVAL OF RADIAL ARTIFACTS

8 CONCLUSIONS AND DISCUSSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: CONVENTION FOR THE FOURIER TRANSFORM, HILBERT TRANSFORM AND CROSS-CORRELATION

APPENDIX B: THEORETICAL REPRESENTATION OF MCBF

APPENDIX C: ARF AND WAVEFIELD SAMPLING

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only