Comparative analysis of formulations for crossed artifact elimination in the frequency–Bessel spectrograms

Abstract

The frequency–Bessel transformation method has significantly advanced the extraction of multimodal surface waves in seismic research. However, the presence of crossed artifacts in frequency–Bessel spectrograms, particularly when stations are regularly distributed, presents a persistent challenge. Various methods have been proposed to mitigate these artifacts, yet their diverse formulations often lead to confusion about their practical application and interrelations. This study aims to demystify these ambiguities by analyzing the existing formulations within a unified framework. We uncover that the apparent discrepancies among these methods primarily originate from the differing conventions across various studies. Consequently, we establish explicit mathematical relationships among these existing formulations. Moreover, we demonstrate that the reliance on numerical Hilbert transform can be avoided by maintaining only the causal component of cross-correlation functions. This approach simplifies the artifact removal process, enhancing the practical utility of frequency–Bessel spectrograms in geophysical analysis.

1 Introduction

Since the introduction of the frequency–Bessel transformation (FJ) method to extract multimodal surface wave dispersion curves by Wang et al. (2019) and Forbriger (2003), it is shortly recognized that crossed artifacts appear in the spectrograms when the stations are regularly distributed. Forbriger (2003) suggested that these artifacts could be reduced to use the Hankel function instead of the Bessel function in the integration process. Following this, subsequent studies by Xi et al. (2021), Luo et al. (2022), Zhou et al. (2023), and Yang et al. (2024) have explored diverse strategies to address this issue.

However, Yang et al. (2024) highlighted that there is a lack of consistency among these methods. For example, Forbriger (2003) and Xi et al. (2021) used in their formulations, while Luo et al. (2022) and Yang et al. (2024) preferred . Additionally, Luo et al. (2022) did not explicitly include the Hilbert transform, which is a requirement in other studies. This disparity raises questions about which method is most effective in practice, whether these methods are fundamentally similar, and how they are related.

In this paper, we show how to extend Forbriger (2003)’s formulation from active to passive sources. We then analyze different approaches, focusing on the works by Forbriger (2003), Xi et al. (2021), Luo et al. (2022), Zhou et al. (2023), and Yang et al. (2024), using a consistent convention, and explore the relationships among these studies. We also make some clarifications on our previous study (Luo et al., 2022) and show how the usage of the numerical Hilbert transform can be avoided.

2 Theory

2.1 Conventions

We use the Hilbert transform in the form (Hahn, 1996)where the symbol P denotes the Cauchy principal value. This definition aligns with current literature (Hahn, 1996; Johnansson, 1999) and is implemented in the MATLAB software and the Python function scipy.signal.hilbert. However, it differs in sign from the definition of Erdelyi (1954) and the Python function scipy.fftpack.hilbert. According to Hu (1989), as cited by Zhou and Chen (2021) and Yang et al. (2024), the Kramers–Kronig relation for an analytic signal is

In this notation, (Erdelyi, 1954, P. 254). Note that there will be a sign difference in Equations 2, 3 if the other definition of the Hilbert transform is used.

2.2 Extend Forbriger's formulation from active to passive source

Because we focus on the component cross-correlation functions (CCFs) in this study, it is shortened as for convenience. The corresponding analytic signal, with as the imaginary part, is denoted as . For the component CCFs , the formulation for the FJ method is as follows: Wang et al. (2019)

In a related study, Forbriger (2003) proposed to use the following equation to mitigate the crossed artifacts:

Note that the original version of is intended for active source seismic records . When applying Equation 5 to ambient noise CCFs, needs to be constructed from the CCFs , which corresponds to the imaginary part of (Sánchez-Sesma and Campillo, 2006; Wang et al., 2019). Using the Kramers–Kronig relation in Equation 2, we obtain

Consequently, for CCFs can be reformulated as

Note that the extension of Forbriger (2003)’s work from active to passive sources intuitively leads to the adoption of .

2.3 Formulation without the numerical Hilbert transform

The FJ method has been extended to multicomponent CCFs (Hu et al., 2020). Correspondingly, the method of Forbriger (2003) to remove crossed artifacts has been adapted for multicomponent CCFs in another slightly different approach by Luo et al. (2022). At this point, we would like to clarify specific aspects of the study by Luo et al. (2022).

Notably, in the approach by Luo et al.’s (2022), there is an implicit assumption that solely the causal part of the CCFs, denoted as , is utilized. Explicitly, . Therefore, the Fourier transform (denoted by ) of will be

Using the Fourier transform of the signum function, and the Fourier transform of a product of functions, , where the symbol represents the convolution and , as well as the definition of the Hilbert transform in Equation 1, there is a relation

From Equation 9, there are relations and . Therefore, . We emphasize that the abovementioned relation offers an alternative approach to computing the Hilbert transform without relying on pure numerical methods. Specifically, it involves first constructing the Fourier transform of the causal part of the CCFs and then extracting the imaginary part of the result. Equation 9 has been validated through numerical tests, as depicted in Figure 1, where is compared with . Specifically, as shown in Figure 1A, we first constructed a causal signal for and then calculated its Fourier transform . The plot juxtaposes in green with in red, demonstrating their concordance and thereby validating Equation 9. Figure 1B depicts a parallel test using CCFs between two virtual stations in a numerical experiment. The phase agreement between the green and red lines is evident. However, a discrepancy in amplitude between the two lines is observed, which is most likely a consequence of edge effects associated with the numerical Hilbert transform.

FIGURE 1

According to Equation 9, the numerical execution of the Hilbert transform can be circumvented in practical scenarios. This is because the Hilbert transform is inherently incorporated during the retention of only the causal portion of the CCFs. Practically, due to the potential asymmetry between the causal and acausal parts of the CCFs, we define the causal part as the average of the positive and negative segments of the CCFs. Additionally, in Luo et al. (2022), only the real part of the integration is used to construct the F-J spectrograms. Therefore, an equivalent expression of the equations used by Luo et al. (2022) would be (see Appendix in Luo et al. (2022)):

Instead of computing Equation 11 directly, Equation 10 indicates that the numerical Hilbert transform can be avoided in the integration if only the casual part of the CCFs is used. Note that in the main text of Luo et al. (2022), in Equations 7, 8, 11, 12 should be replaced by . It is also noteworthy that although the usage of Re in Equation 10 is not mentioned in Luo et al. (2022), it is required because the corresponding imaginary part of the integrand (denoted by ) fails to correctly constrain the dispersion curves (see details in the Discussion section).

2.4 Comparison among different formulations

In the following section, we summarize other formulations for the improved FJ method with crossed artifacts removed.

Xi et al. (2021) proposed to usewhere is the analytic signal defined in Equation 6 and is the complex conjugate of .

Zhou et al. (2023) proposed to usewhere should be regarded as the analytic signal whose real part is the CCFs . Therefore, and according to Equation 3. Note that , and Equation 13 seems to be a correct version of Zhou and Chen’s (2021)’s original formulation.

The methods of Xi et al. (2021) and Zhou et al. (2023) have been implemented by Li et al. (2021) for ZZ component CCFs and by Zhang et al. (2022) for multicomponent CCFs. It is important to note, however, that Li et al. (2021) and Zhang et al. (2022) utilized the Python function scipy. fftpack.hilbert in their implementation. This usage has introduced an apparent sign difference when compared to the theoretical analysis presented in this paper.

In a recent study, Yang et al. (2024) proposed the frequency–Hankel transform

Combining with the properties of the Hankel function, , , and , formulations presented in Equations 7, 10, 12–14 can be expanded in terms of Bessel functions asand

According to Equations 15–19, there relation is given as

Denoting the corresponding imaginary part of the integrand in Equation 10 as , the relation is given as

Note that Equation 19 has already been correctly provided by Yang et al. (2024) (Equation 21). However, the construction of in Yang et al. (2024) seems to be problematic (it should be instead of on page KS72 in Yang et al. (2024) in their notation according to our Equation 6). The problematic derivation results in their flawed arguments regarding the relationships among different methods.

3 Discussion

As elucidated in Equations 20, 21, the formulations proposed by Forbriger (2003), Xi et al. (2021), Luo et al. (2022), Zhou et al. (2023), and Yang et al. (2024) are actually correlated with each other. Notably, the constructions of , , , and all require the numerical computation of the Hilbert transform, as implemented by Li et al. (2021) and Zhang et al. (2022). We would like to emphasize that the numerical Hilbert transform can be avoided in the implementation of by merely maintaining the causal portion of the CCFs, or pragmatically, by designating the average of the positive and negative parts of the CCFs as the causal part. To substantiate the efficacy of , validate the relationship in Equation 20, and study the effects of the real part of and , we present a numerical example. In practice, Li et al. (2021) and Zhang et al. (2022) implemented and using Equations 16, 18, instead of their original formulations in Equations 12, 13. We, therefore, choose to compare the results of and using Equations 10, 16 (Figures 2B, C). We also plot to study the effects of the real parts of and (Figure 2D). In the numerical example, the seismic arrays are located linearly, as shown in Figure 2B in Luo et al. (2022), and the velocity model contains a low-velocity layer as shown in Figure 2C in Luo et al. (2022). The synthetic waveforms were computed using the Computer Programs in Seismology package (Herrmann, 2013). Subsequently, the vertical–vertical component CCFs were computed. Based on this, the FJ spectrograms of the synthetic data, in accordance with Equations 4, 10, 16, were calculated and illustrated in Figures 2A–C, respectively. In this example, it is clear the crossed artifacts are effectively removed in both Figures 2B, C. Notably, the outputs of and are almost identical to each other, with the exception that the normalization factor of at each frequency is four times that of , as delineated in Equation 20. On the other hand, the peak locations of the spectrogram, corresponding to Re and Re, do not coincide with the correct dispersion curves (Figure 2D). Therefore, as discussed by Zhou and Chen (2021), the inclusion of Re reduces the resolution of FJ spectrograms.

FIGURE 2

To further demonstrate the correlations among these methods, we conducted a parallel test using field data from a USArray, as shown in Supplementary Figure S1. The data were previously used in Luo et al. (2022) (see their Figure 9). This test also confirms that both and are equally effective in removing the crossed artifacts. Additionally, the real part of should not be included when generating high-resolution FJ spectrograms.

4 Conclusion

In this study, we have conducted analysis of prior work addressing the elimination of crossed artifacts in the frequency–Bessel spectrograms, placing these studies within a unified framework. This approach has allowed us to uncover the underlying relationships between various methodologies proposed in the literature. We also have shown how to extend the formulation proposed by Forbriger (2003) from active to passive sources. Furthermore, we have proposed a refined approach that circumvents the need for the numerical Hilbert transform, as has been implicitly suggested in a previous study (Luo et al., 2022).

References

• 1

  ErdelyiA. (1954). Tables of integral transforms, vol. 2, bateman manuscript project. New York, New York: McGraw-Hill.

  • Google Scholar

• 2

  ForbrigerT. (2003). Inversion of shallow-seismic wavefields: I. wavefield transformation. Geophys. J. Int.153, 719–734. 10.1046/j.1365-246x.2003.01929.x

  • CrossRef
  • Google Scholar

• 3

  HahnS. L. (1996). Hilbert transforms in signal processing. Boston, London: Artech House.

  • Google Scholar

• 4

  HerrmannR. B. (2013). Computer programs in seismology: an evolving tool for instruction and research. Seismol. Res. Lett.84, 1081–1088. 10.1785/0220110096

  • CrossRef
  • Google Scholar

• 5

  HuB. Y. K. (1989). Kramers–Kronig in two lines. Am. J. Phys.57, 821. 10.1119/1.15901

  • CrossRef
  • Google Scholar

• 6

  HuS.LuoS.YaoH. (2020). The frequency-Bessel spectrograms of multicomponent cross-correlation functions from seismic ambient noise. J. Geophys. Res. Solid Earth125, e2020JB019630. 10.1029/2020JB019630

  • CrossRef
  • Google Scholar

• 7

  JohnanssonM. (1999). The Hilbert transform. Master’s thesis. Suecia: Växjö University.

  • Google Scholar

• 8

  LiZ.ZhouJ.WuG.WangJ.ZhangG.DongS.et al (2021). CC-FJpy: a python package for extracting overtone surface-wave dispersion from seismic ambient-noise cross correlation. Seismol. Res. Lett.92, 3179–3186. 10.1785/0220210042

  • CrossRef
  • Google Scholar

• 9

  LuoS.HuS.ZhouG.YaoH. (2022). Improvement of frequency-Bessel phase-velocity spectra of multicomponent cross-correlation functions from seismic ambient noise. Bull. Seismol. Soc. Am.112, 2257–2279. 10.1785/0120220027

  • CrossRef
  • Google Scholar

• 10

  Sánchez-SesmaF. J.CampilloM. (2006). Retrieval of the Green’s function from cross correlation: the canonical elastic problem. Bull. Seismol. Soc. Am.96, 1182–1191. 10.1785/0120050181

  • CrossRef
  • Google Scholar

• 11

  WangJ.WuG.ChenX. (2019). Frequency-Bessel transform method for effective imaging of higher-mode Rayleigh dispersion curves from ambient seismic noise data. J. Geophys. Res. Solid Earth124, 3708–3723. 10.1029/2018JB016595

  • CrossRef
  • Google Scholar

• 12

  XiC.XiaJ.MiB.DaiT.LiuY.NingL. (2021). Modified frequency-Bessel transform method for dispersion imaging of Rayleigh waves from ambient seismic noise. Geophys. J. Int.225, 1271–1280. 10.1093/gji/ggab008

  • CrossRef
  • Google Scholar

• 13

  YangZ.SunY.-C.ZhangD.HanP.ChenX. (2024). A frequency-Hankel transform method to extract multimodal Rayleigh wave dispersion spectra from active and passive source surface wave data. Geophysics89, KS69–KS81. 10.1190/geo2023-0189.1

  • CrossRef
  • Google Scholar

• 14

  ZhangG.LiuQ.ChenX. (2022). Enhancing the frequency-Bessel spectrogram of ambient noise cross-correlation functions. Bull. Seismol. Soc. Am.113, 361–377. 10.1785/0120220124

  • CrossRef
  • Google Scholar

• 15

  ZhouJ.ChenX. (2021). Removal of crossed artifacts from multimodal dispersion curves with modified frequency-Bessel method. Bull. Seismol. Soc. Am.112, 143–152. 10.1785/0120210012

  • CrossRef
  • Google Scholar

• 16

  ZhouJ.LiZ.ChenX. (2023). Extending the frequency band of surface wave dispersion curves by combining ambient noise and earthquake data and self-adaptive normalization. J. Geophys. Res. Solid Earth128, e2022JB026040. 10.1029/2022JB026040

  • CrossRef
  • Google Scholar