**1. Introduction**

Passive seismic methods do not require the excitation of explosives, vibroseis, etc.; they simply require the placement of node geophones according to the active source acquisition line and point spacing, and the data can be collected by continuously receiving ambient noises. Then, by using seismic interferometry (SI), the responses of virtual sources can be retrieved through the correlated responses on the two receivers [1–4]. The term "interferometry" generally refers to the study of interference between signals in order to obtain information through the differences between them. In the field of SI, it is used to study the interference of seismic-related signals. Its basic operation steps are very simple, mainly divided into two steps: the first step is cross-correlation calculation, which can be understood as detecting the travel time difference of waves recorded between two receivers; the second step is stacking—that is, integrating all actual sources.

Although the theory of SI is not limited to either body or surface waves [5,6], extracting surface waves from the cross-correlation of ambient noise is quite robust and has become a mature technology [7]. In contrast, it is much more difficult to extract body waves [8,9]. However, some studies have also shown its possibility at various scales [10–13]. Even so, in most SI studies, retrieving surface waves is still one of the major applications. Surface waves with different periods can be used to investigate the earth structure at different depths [14]. Generally speaking, most studies focus on frequencies lower than 1 Hz to provide images of the structure of the crust and upper mantle [15,16]. However, some studies have proven that structures tens to hundreds of meters underground can be imaged with frequencies higher than 1 Hz, such as by using traffic noise to detect urban underground space [17,18], which has attracted the attention of more and more engineering seismologists.

**Citation:** Fang, J.; Liu, G.; Liu, Y. Application of an Automatic Noise or Signal Removal Algorithm Based on Synchrosqueezed Continuous Wavelet Transform of Passive Surface Wave Imaging: A Case Study in Sichuan, China. *Appl. Sci.* **2021**, *11*, 11718. https://doi.org/10.3390/ app112411718

Academic Editor: Amerigo Capria

Received: 6 October 2021 Accepted: 7 December 2021 Published: 9 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Green's function reconstruction of diffuse wavefield (i.e., waves arrive from all angles with equal strength) observations in an open configuration is the basis for SI with passive data [19]. However, it is not realistic to obtain a pure diffuse wavefield from observations. Non-ideal source distributions such as limited azimuthal distribution or near-source effects [20] can significantly affect the reliability and interpretation of observations due to the violation of the stationarity assumptions required by passive survey methods [21]. In order to improve the accuracy of the dispersion measurement, some studies have begun to consider how to attenuate the influence of noise sources with prominent directivity. Aki [5] proposed a classical spatial autocorrelation method which used the mathematical transformation of a symmetric receiving array (e.g., circle) to eliminate the azimuth factor of the noise source. Park et al. [22] introduced a novel strategy for imaging dispersion of passive surface waves with an active scheme based on phase-shift measurement, called roadside passive multichannel analysis of surface waves. Scanning processes were implemented along potential azimuthal directions of noise. Feuvre et al. [23] proposed a cross-correlation-based beamforming average dispersion analysis method which effectively restrained the directionality of noise and reduced the aliasing effect. Cheng et al. [24] proposed a multichannel analysis of the passive surface (MAPS) wave method based on long noise sequence cross-correlation to handle the azimuthal effects for directional noise sources. Distinct from the above methods, Mousavi et al. [25] proposed a new and fast algorithm for accurate noise removal/signal removal based on higher-order statistics (HOS), general cross-validation (GCV), and wavelet hard thresholding (WHT) in synchrosqueezed domains and tested the performance of the proposed algorithm with synthetic and real seismic data. In addition, he indicated that it can also be an effective procedure in ambient noise studies.

In this paper, we applied the noise or signal removal algorithm based on SS-CWT proposed by Mousavi et al. [25] to a 1-h passive seismic dataset acquired in Sichuan, China, and verified the effectiveness of this method in enhancing passive surface waves from short time series. First, we used this method to preprocess the 1-hour passive seismic data in Sichuan to remove the prominent noise events and used cross-correlation to produce virtual shot gathers from the preprocessed noise records. Then, after cross-correlation, the method was used as a post-processing step to denoise the virtual shot gathers, improving the signalto-noise ratio (S/N) of the virtual shot gathers and further improving the reconstructed surface waves. Next, after preprocessing and post-processing, virtual shot gathers were used for dispersion analysis with an active scheme based on the phase-shift measurement, and inversion was conducted to obtain the underground shear wave velocity profile in this area. Finally, the applicability and shortcomings of this method are discussed.

### **2. Methodology**

The methodology used in this paper was proposed by Mousavi et al. [25]. We summarized and sorted out their methods, mainly including the following three processing steps.

### *2.1. Preprocessing*

For the passive seismic data *y* collected in the field, first, we carried out CWT and calculated its coefficient *Wy*, which represents the finite energy of the data in a concentrated time–frequency picture. CWT is a multiresolution transform method given by [26], as shown in Equation (1):

$$\mathcal{W}\_{\mathcal{Y}}(a,\tau) = \langle \mathcal{y}, \psi\_{a,\tau} \rangle = \int\_{-\infty}^{+\infty} y(t) a^{-\frac{1}{2}} \psi^\* \left( \frac{t-\tau}{a} \right) dt,\tag{1}$$

where *a* and *τ* denote the scale and time shift, respectively; ∗ is the complex conjugate; *y*, *ψ* denotes the inner product; and *ψ* is the mother wavelet.

Then, we used HOS and the kurtosis criterion to detect the scale of and remove the pure Gaussian noise correlation coefficient from the time–frequency analysis (TFR), leaving the combination of noise and signal. The *Kurt* of *N* observation coefficients *Wy* can be calculated by Equation (2):

$$Kurt\_y = \frac{\sum\_{n=1}^{N} \left(W\_{y\_n} - \mu\_{W\_y}\right)^4}{N\sigma\_{W\_y}^4} - 3,\tag{2}$$

where *σWy* and *μWy* are the estimated standard deviation and mean value of wavelet coefficient *Wy*, respectively, and *N* denotes the sample sequence. HOS measure the degree of Gaussianity, so they have been used as detectors of non-Gaussian signals in Gaussian noise [27]. Equation (3) defines the HOS criterion, which is retained for the Gaussianity measure for distinguishing Gaussian distributions from non-Gaussian distributions:

$$|Kurt\_y| \le \frac{\sqrt{24/N}}{\sqrt{1-a}}\text{.}\tag{3}$$

where *α* is the level of confidence. Ravier et al. [28] estimated that the optimal value of *α* is 90%.

### *2.2. Thresholding*

After obtaining the processed coefficient *Wy* through the preprocessing steps, the synchrosqueezed transform was performed according to Equation (4) to obtain the synchrosqueezed-CWT (SS-CWT) coefficients *Ty*.

$$T\_{\mathcal{Y}}(\omega\_{\ell\prime}\tau) = \Delta\omega^{-1} \sum\_{\substack{a\_k \mid \omega\left(a\_k,\tau\right) - \omega\_\ell\left| \leq \frac{\Delta\omega}{2}\right.}} W\_{\mathcal{Y}}(a\_k,\tau)a\_k^{-\frac{3}{2}}\Delta a\_{k\prime} \tag{4}$$

where *ω* is the discrete frequency and *ω* represents the th discrete frequency, *ak* represents the *k*th scale, *τ* is the time shift, and Δ represents the increment symbol, where Δa = *ak* − *ak*−1.

Then, we used the hard threshold scheme (Equation (5)) proposed by Donoho et al. [29,30] to provide a threshold for the SS-CWT coefficients *Ty*.

$$\eta\_{\lambda}^{h}(\mathcal{W}\_{\mathcal{Y}}\lambda) = \left\{ \begin{array}{c} \mathcal{W}\_{\mathcal{Y}} \ \text{ } \prime \left| \mathcal{W}\_{\mathcal{Y}} \right| > \lambda \\ 0 \end{array} \right. \\ \left. \begin{array}{c} \left| \mathcal{W}\_{\mathcal{Y}} \right| > \lambda \\ \left| \mathcal{W}\_{\mathcal{Y}} \right| < \lambda \end{array} \right. \tag{5}$$

where λ is the selected appropriate threshold; η is the selected threshold rule, which can be divided into a hard threshold and a soft threshold; η*<sup>h</sup>* <sup>λ</sup> represents the hard threshold rule; and *Wy* represents the wavelet coefficients of observation *y*.

The optimal threshold λ was automatically determined by the GCV method, which was proposed by Nason et al. [31]. The GCV function is defined as follows:

$$\text{GCV}(\lambda) = \frac{1}{N} \frac{||T\_{\mathcal{Y}} - \bar{T}\_{\lambda}||^2}{||\frac{N\_0}{N}||^2},\tag{6}$$

where *T* -<sup>λ</sup> is the threshold coefficient using the threshold λ, and *N*<sup>0</sup> is the number of coefficients using the threshold λ to return to zero. This function simulated the error between the estimated signal and the real signal, so its minimum value was used to select an optimal threshold. By selecting thresholds for the main components of the signal, the initial estimation of the signal was obtained using Equation (7):

$$y\_k(t) = 2C\_\psi^{-1} Re\left(\sum\_{l \in L\_k(t)} T\_\mathcal{Y}(\omega\_{\ell'}, t)\right),\tag{7}$$

where *Lk*(*t*) denotes a small frequency band around the ridge of the *k*th component in the SS-CWT, the symbol *Re* represents the real part, and the constant *C<sup>ψ</sup>* is given by Thakur et al. [32]:

$$\mathbb{C}\_{\Psi} = \int\_{0}^{\infty} \xi^{-1} \hat{\Psi}^\*(\xi) d\xi,\tag{8}$$

where ∗ represents the complex conjugate; *ψ* is the mother wavelet, which is a square integrateable function; *ξ* is the angular frequency; and *ψ*ˆ is the Fourier transform of *ψ*.
