Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1505367 10.3389/feart.2025.1505367 Earth Science Mini Review Comparative analysis of formulations for crossed artifact elimination in the frequency–Bessel spectrograms Jiang et al. 10.3389/feart.2025.1505367 Jiang Xiaohuan 1 2 Hu Shaoqian 2 * Luo Song 3 1 School of Civil Engineering and Architecture, Wuhan Polytechnic University, Wuhan, China 2 Hubei Subsurface Multi-Scale Imaging Key Laboratory, School of Geophysics and Geomatics, China University of Geosciences, Wuhan, China 3 School of Geological Engineering and Geomatics, Chang’an University, Xi’an, China

Edited by: Elsa Leticia Flores Márquez, National Autonomous University of Mexico, Mexico

Reviewed by: Kai Zhang, Huairou Lab, China

 *Correspondence: Shaoqian Hu, hushaoqian@cug.edu.cn 04 03 2025 2025 13 1505367 02 10 2024 12 02 2025 Copyright © 2025 Jiang, Hu and Luo. 2025 Jiang, Hu and Luo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

 multimodal surface wave frequency–Bessel transformation crossed artifact elimination Hilbert transform comparative analysis section-at-acceptance Solid Earth Geophysics 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 H 0 ( 2 ) in their formulations, while Luo et al. (2022) and Yang et al. (2024) preferred H 0 ( 1 ) . 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) H f ω = 1 π P ∫ − ∞ ∞ f ω ′ ω − ω ′ d ω ′ , 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 χ ( ω ) = χ 1 ( ω ) + i χ 2 ( ω ) is χ 1 ω = − H χ 2 ω , χ 2 ω = H χ 1 ω .

In this notation, H ( sgn ( w ) J 0 ( w ) ) = Y 0 ( | w | ) (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 Z Z component cross-correlation functions (CCFs) C Z Z ( ω , r ) in this study, it is shortened as C ( ω , r ) for convenience. The corresponding analytic signal, with C ( ω , r ) as the imaginary part, is denoted as G ( ω , r ) . For the Z Z component CCFs C ( ω , r ) , the formulation for the FJ method is as follows: Wang et al. (2019) I Wang ω , k = ∫ 0 ∞ C ω , r J 0 k r r d r .

In a related study, Forbriger (2003) proposed to use the following equation to mitigate the crossed artifacts: I Forbriger 0 ω , k = ∫ 0 ∞ G ω , r H 0 2 k r r d r .

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

Consequently, I Forbriger for CCFs can be reformulated as I Forbriger ω , k = ∫ 0 ∞ − H C ω , r + i C ω , r H 0 2 k r d r .

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

 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 c ̄ ( t ) , is utilized. Explicitly, c ̄ ( t ) = 1 2 c ( t ) + sgn ( t ) c ( t ) . Therefore, the Fourier transform (denoted by F ) of 2 c ̄ ( t ) will be 2 C ̄ ω = F c t + sgn t c t = C ω − i F i sgn t c t ω .

Using the Fourier transform of the signum function, F ( i sgn ( t ) ) = 2 ω and the Fourier transform of a product of functions, F [ g ( t ) f ( t ) ] = 1 2 π g ̂ ∗ f ̂ ( ω ) , where the symbol ∗ represents the convolution and F ( f ) = f ̂ ( ω ) , as well as the definition of the Hilbert transform in Equation 1, there is a relation 2 C ̄ ω = C ω − i 1 π 1 ω ∗ C ω = C ω − i H C ω .

From Equation 9, there are relations Re { C ̄ ( ω ) } = 1 2 C ( ω ) and Im { C ̄ ( ω ) } = − 1 2 H [ C ] ( ω ) . Therefore, H [ C ] ( ω ) = − 2 Im { C ̄ ( ω ) } . 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 ( c ̄ ( t ) ) and then extracting the imaginary part of the result. Equation 9 has been validated through numerical tests, as depicted in Figure 1, where Im { C ̄ ( ω ) } is compared with − 1 2 H [ C ] ( ω ) . Specifically, as shown in Figure 1A, we first constructed a causal signal c ̄ ( t ) = 1 2 e − t 2 for t > 0 and then calculated its Fourier transform C ̄ ( ω ) . The plot juxtaposes Im { C ̄ ( ω ) } in green with − 1 2 H [ C ] ( ω ) 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.

Validation of Equation 9 by comparison of different methods to construct the imaginary part of C ̄ ( ω , r ) : (A) with the input signal described as c ̄ ( t ) = e − t 2 / 2 for t > 0 , (B) with the input signal as a sample CCF between two virtual stations in a numerical test. Blue lines are the real part of C ̄ ( ω , r ) , which equals 1 2 C ( ω , r ) . Green lines are the imaginary part of C ̄ ( ω , r ) constructed by the Fourier transform of the causal part of CCFs ( Im { C ̄ ( ω ) } ) . Red lines are the imaginary part of C ̄ ( ω , r ) constructed by the numerical Hilbert transform − 1 2 H [ C ] ( ω ) .

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)): I Luo ω , k = Re ∫ 0 ∞ C ̄ ω , r H 0 1 k r r d r , = 1 2 Re ∫ 0 ∞ C ω , r − i H C ω , r H 0 1 k r r d r .

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 C ̄ ( ω , r ) is used. Note that in the main text of Luo et al. (2022), H 0 ( 2 ) in Equations 7, 8, 11, 12 should be replaced by H 0 ( 1 ) . 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 I * Luo ) 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 use I Xi ω , k = − ∫ 0 ∞ i G ω , r H 0 2 k r + i G ω , r * H 0 1 k r r d r , where G is the analytic signal defined in Equation 6 and ( i G ) * is the complex conjugate of i G .

Zhou et al. (2023) proposed to use I Zhou ω , r = ∫ 0 ∞ Re C ̃ r , ω J 0 k r + Im C ̃ r , ω Y 0 k r r d r , where C ̃ ( ω , r ) should be regarded as the analytic signal whose real part is the CCFs C ( ω , r ) . Therefore, Re C ̃ ( ω , r ) = C ( ω , r ) and Im C ̃ ( ω , r ) = H [ C ] ( ω , r ) according to Equation 3. Note that C ̃ ( ω , r ) ≠ G ( ω , r ) , 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 I Yang ω , k = ∫ 0 ∞ H C ω , r + i C ω , r H 0 1 k r r d r .

Combining with the properties of the Hankel function, H 0 ( 1 ) ( x ) = J 0 ( x ) + i Y 0 ( x ) , H 0 ( 2 ) ( x ) = J 0 ( x ) − i Y 0 ( x ) , and ( H 0 1 ( x ) ) * = H 0 ( 2 ) ( x ) , formulations presented in Equations 7, 10, 12–14 can be expanded in terms of Bessel functions as I Forbriger = ∫ 0 ∞ − H C ω , r J 0 k r + C ω , r Y 0 k r r d r + i ∫ 0 ∞ H C ω , r Y 0 k r + C ω , r J 0 k r r d r , I Xi = 2 ∫ 0 ∞ H C ω , r Y 0 k r + C ω , r J 0 k r r d r , I Luo = 1 2 ∫ 0 ∞ H C ω , r Y 0 k r + C ω , r J 0 k r r d r , I Zhou = ∫ 0 ∞ H C ω , r Y 0 k r + C ω , r J 0 k r r d r , and I Yang = ∫ 0 ∞ H C ω , r J 0 k r − C ω , r Y 0 k r r d r + i ∫ 0 ∞ H C ω , r Y 0 k r + C ω , r J 0 k r r d r .

According to Equations 15–19, there relation is given as Im I Forbriger = 1 2 I Xi = 2 I Luo = I Zhou = Im I Yang .

Denoting the corresponding imaginary part of the integrand in Equation 10 as I * Luo , the relation is given as Re I Forbriger = 2 I * Luo = − Re I Yang .

Note that Equation 19 has already been correctly provided by Yang et al. (2024) (Equation 21). However, the construction of G ( ω , r ) in Yang et al. (2024) seems to be problematic (it should be − a + b i instead of a + b i 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 I Forbriger , I Xi , I Zhou , and I Yang 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 I Luo 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 I Luo , validate the relationship in Equation 20, and study the effects of the real part of I Forbriger and I Yang , we present a numerical example. In practice, Li et al. (2021) and Zhang et al. (2022) implemented I Xi and I Zhou using Equations 16, 18, instead of their original formulations in Equations 12, 13. We, therefore, choose to compare the results of I Luo and I Xi using Equations 10, 16 (Figures 2B, C). We also plot I * Luo to study the effects of the real parts of I Forbriger and I Yang (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 I Luo and I Xi are almost identical to each other, with the exception that the normalization factor of I Xi at each frequency is four times that of I Luo , as delineated in Equation 20. On the other hand, the peak locations of the I * Luo spectrogram, corresponding to Re { I Forbriger } and Re { I Yang } , do not coincide with the correct dispersion curves (Figure 2D). Therefore, as discussed by Zhou and Chen (2021), the inclusion of Re { I Forbriger } reduces the resolution of FJ spectrograms.

Frequency–Bessel spectrograms of synthetic data computed by the FJ method and three modified versions by Luo et al. (2022) and Xi et al. (2021). (A) FJ spectrogram by integrating with the Bessel function J 0 as I Wang in Equation 4. (B) FJ spectrogram by I Luo as expressed in Equation 10. (C) FJ spectrogram by I Xi as expressed in Equation 16. (D) FJ spectrogram by I Luo * , which is the imaginary part of the corresponding integrand in Equation 10. The numerical Hilbert transform was performed to construct I Xi . The theoretical dispersion curves are plotted as red crosses. The normalization factors at each frequency are also plotted at the top of the FJ spectrograms (divided by 1 0 4 for visualization).

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 I Luo and I Xi are equally effective in removing the crossed artifacts. Additionally, the real part of I Forbriger 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).

 Author contributions

XJ: formal analysis, investigation, validation, visualization, and writing–original draft. SH: conceptualization, formal analysis, methodology, validation, and writing–original draft, Writing–review and editing. SL: Software, Writing–review and editing.

 Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This study was supported by the National Natural Science Foundation of China (42488301, 41904047), Hubei Subsurface Multi-scale Imaging Key Laboratory (China University of Geosciences) (SMIL-2023-8), and the Research and Innovation Initiatives of WHPU.

The authors would like to thank the reviewers for their comments that helped them improve the manuscript.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2025.1505367/full#supplementary-material

References Erdelyi A. (1954). Tables of integral transforms, vol. 2, bateman manuscript project. New York, New York: McGraw-Hill. Forbriger T. (2003). Inversion of shallow-seismic wavefields: I. wavefield transformation. Geophys. J. Int. 153, 719–734. 10.1046/j.1365-246x.2003.01929.x Hahn S. L. (1996). Hilbert transforms in signal processing. Boston, London: Artech House. Herrmann R. B. (2013). Computer programs in seismology: an evolving tool for instruction and research. Seismol. Res. Lett. 84, 1081–1088. 10.1785/0220110096 Hu B. Y. K. (1989). Kramers–Kronig in two lines. Am. J. Phys. 57, 821. 10.1119/1.15901 Hu S. Luo S. Yao H. (2020). The frequency-Bessel spectrograms of multicomponent cross-correlation functions from seismic ambient noise. J. Geophys. Res. Solid Earth 125, e2020JB019630. 10.1029/2020JB019630 Johnansson M. (1999). The Hilbert transform. Master’s thesis. Suecia: Växjö University. Li Z. Zhou J. Wu G. Wang J. Zhang G. Dong S. (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 Luo S. Hu S. Zhou G. Yao H. (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 Sánchez-Sesma F. J. Campillo M. (2006). Retrieval of the Green’s function from cross correlation: the canonical elastic problem. Bull. Seismol. Soc. Am. 96, 1182–1191. 10.1785/0120050181 Wang J. Wu G. Chen X. (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. 10.1029/2018JB016595 Xi C. Xia J. Mi B. Dai T. Liu Y. Ning L. (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 Yang Z. Sun Y.-C. Zhang D. Han P. Chen X. (2024). A frequency-Hankel transform method to extract multimodal Rayleigh wave dispersion spectra from active and passive source surface wave data. Geophysics 89, KS69–KS81. 10.1190/geo2023-0189.1 Zhang G. Liu Q. Chen X. (2022). Enhancing the frequency-Bessel spectrogram of ambient noise cross-correlation functions. Bull. Seismol. Soc. Am. 113, 361–377. 10.1785/0120220124 Zhou J. Chen X. (2021). Removal of crossed artifacts from multimodal dispersion curves with modified frequency-Bessel method. Bull. Seismol. Soc. Am. 112, 143–152. 10.1785/0120210012 Zhou J. Li Z. Chen X. (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 Earth 128, e2022JB026040. 10.1029/2022JB026040