EEMD and Multiscale PCA-Based Signal Denoising Method and Its Application to Seismic P-Phase Arrival Picking

Abstract

Signal denoising is one of the most important issues in signal processing, and various techniques have been proposed to address this issue. A combined method involving wavelet decomposition and multiscale principal component analysis (MSPCA) has been proposed and exhibits a strong signal denoising performance. This technique takes advantage of several signals that have similar noises to conduct denoising; however, noises are usually quite different between signals, and wavelet decomposition has limited adaptive decomposition abilities for complex signals. To address this issue, we propose a signal denoising method based on ensemble empirical mode decomposition (EEMD) and MSPCA. The proposed method can conduct MSPCA-based denoising for a single signal compared with the former MSPCA-based denoising methods. The main steps of the proposed denoising method are as follows: First, EEMD is used for adaptive decomposition of a signal, and the variance contribution rate is selected to remove components with high-frequency noises. Subsequently, the Hankel matrix is constructed on each component to obtain a higher order matrix, and the main score and load vectors of the PCA are adopted to denoise the Hankel matrix. Next, the PCA-denoised component is denoised using soft thresholding. Finally, the stacking of PCA- and soft thresholding-denoised components is treated as the final denoised signal. Synthetic tests demonstrate that the EEMD-MSPCA-based method can provide good signal denoising results and is superior to the low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, EEMD-Hankel–SVD, and wavelet-MSPCA-based denoising methods. Moreover, the proposed method in combination with the AIC picking method shows good prospects for processing microseismic waves.

1. Introduction

Signals from various practical applications are usually corrupted by unwanted noise, such as background noise, electronic noise, and instrumental noise. Consequently, to improve signal detection, classification, and phase arrival picking, it is necessary to remove unwanted noise. Currently, signal denoising is widely used in the fields of medical signals [1], speech signals [2], mechanical diagnosis [3], structural health monitoring [4,5,6], and seismic monitoring networks [7,8,9]; however, this study focused specifically on microseismic (MS) signal denoising.

Spectral filtering is frequently applied to denoise a signal, and a Butterworth-based filter (low-pass filter, high-pass filter, and band-pass filter) is usually adopted [10], which is simple but may filter out a frequency band that contains useful information when the noise and signal share the same frequency band. In addition, the adopted frequency band may vary for different signals; in other words, it cannot denoise a signal adaptively. Therefore, more effective filtering methods have been proposed for signal denoising, especially based on multiscale features (e.g., wavelet analysis and empirical mode decomposition).

For wavelet analysis-based signal denoising, Beenamol et al. [7] used Shannon and Tsallis entropies to select important wavelet coefficients to denoise a signal; however, noise may still exist in the selected wavelet coefficients. Thus, Tang et al. [8], Mousavi et al. [9], and Chen et al. [11] denoised MS signals based on the combination of wavelet transformation and adaptive thresholding, synchro squeezing wavelet transformation and custom thresholding of single-channel data, and empirical wavelet transformation and automatic noise attenuation thresholding algorithms. However, redundant information in the low-frequency band is not removed by filtering out the high-frequency components and thresholding. Patil [12] combined wavelet transformation and singular value decomposition (SVD) for noise reduction, whereas Bakshi [13] decomposed signals by wavelet transformation and applied principal component analysis (PCA) to denoise the wavelet coefficients of multiple signals containing similar noises and obtained good denoising results. SVD and PCA can efficiently separate useful signals from noise. Wavelet analysis decomposes signals from high- to low-frequency components through wavelet function expansion and translation, which has limited adaptive decomposition abilities for complex signals.

Empirical mode decomposition (EMD), proposed by Huang et al. [14], has better adaptability than wavelet decomposition, and EMD and its improvement have been widely used in signal processing. Mao et al. [15], Dao et al. [16], and Li et al. [17] proposed signal denoising methods based on EMD and thresholding; however, some redundant information could not be denoised by thresholding. Qiao et al. [18] adopted EMD with multiscale eigenvalues of SVD for denoising, which could further reduce noise. However, similar components may exist in different scales for EMD, making it difficult to separate useful components and noise. To overcome this problem, Wu and Huang [19] proposed ensemble empirical mode decomposition (EEMD). On this basis, Jia et al. [20], Han et al. [21], and Li et al. [17] took advantage of permutation entropy, adaptive thresholding, and soft thresholding to denoise signals, attaining a better denoising performance than that based on EMD. Moreover, variational mode decomposition (VMD) has been applied for signal denoising [22], which successfully avoids mode mixing and is robust for both noise and sampling, unlike EMD.

Subsequently, SVD was combined with the Hankel matrix to denoise a signal [23], where the Hankel matrix can expand a signal into a high-dimensional matrix and is very beneficial for feature extraction. Recently, coda wave interferometry (CWI) has been successfully applied to monitor stress distribution and damage detection in engineering [4,5,6], where wavelet transformation [5], normal band-pass filter, narrow band-pass filter, and double-notch filter [6] were used to extract useful signal features. These techniques provided a good research direction for noise reduction and damage detection. In addition, deep neural networks have been proposed for signal denoising; for example, DeepDenoiser [24] achieves impressive denoising of seismic signals, even when the signal and noise share a common frequency band. However, perfect separation of signal and noise in the time–frequency domain requires recovering two complex spectra for the signal and noise from the complex spectra of the noisy signal.

In this study, we focused on the PCA denoising concept proposed by Bakshi [13], which can enhance useful signal and noise separation, and the highly adaptive EEMD, Hankel matrix, and multiscale PCA (MSPCA) are combined to remove noise from a single signal. The main innovations and contributions of this research are as follows: (1) The EEMD was introduced into MSPCA denoising to achieve improved signal decomposition adaptability compared to wavelet decomposition. (2) The Hankel matrix was used to transform a single EEMD component into high-dimensional data that are suitable for conducting PCA denoising on only one signal. (3) Soft thresholding was selected for further denoising of the PCA-denoised component. (4) The denoising results are superior to those of the low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, EEMD-Hankel–SVD, and wavelet-MSPCA-based denoising methods. The remainder of this paper is organised as follows. Section 2 introduces the EEMD, Hankel matrix, and MSPCA-based denoising methods. The synthetic and application tests are described in Section 3 and Section 4, respectively. Section 5 compares the proposed method with other denoising methods. The research results and prospects are presented in Section 6.

2. Theoretical Fundamentals

2.1. Basis of EEMD

Wu and Huang [19] proposed a noise-assisted EEMD based on EMD by adding white noise to the signals. The steps of the EEMD algorithm are as follows.

(1) By adding a group of white Gaussian noise ω_j(t) into a noisy signal x(t), a new signal $x_j^{\prime }\left(t\right)$ is obtained as follows:

$$x_j^{\prime }\left(t\right)=x\left(t\right)+{\omega }_j\left(t\right)$$

(1)

where t is the signal time, j is the jth time of adding white Gaussian noise, and j = 1, 2,…, M, where M is the number of times white Gaussian noise is added.

(2) Through EMD of $x_{{}_j}^{\prime }\left(t\right)$, n intrinsic mode functions (IMFs), $c_{j,i}$, and one residual component $c_{j,n+1}$ are obtained as follows:

$$x_{{}_j}^{\prime }(t)={\displaystyle \sum _{i=1}^n{c}_{j,i}}(t)+c_{j,n+1}$$

(2)

where i is the IMF identity, i = 1, 2,…, n, and n is the number of IMFs.

(3) By adding different white Gaussian noises [ω_j(t)] into the noisy signal x(t) and repeating steps (1) and (2), the ith IMF component $c_{j,i}(t)$ and the residual component $c_{j,n+1}$ through the jth decomposition are obtained.

(4) Based on the zero-mean principle of the spectra of white Gaussian noise, the effects of white Gaussian noise as the reference structure of the time-domain distribution are eliminated. The average IMF component c_i(t) and the residual component $c_{j,n+1}$ corresponding to noisy signals can be expressed as follows:

$$\{\begin{array}{l}{c}_i\left(t\right)=\frac{1}{M}{\displaystyle \sum _{j=1}^M{c}_{j,i}}(t)\\ c_{n+1}\left(t\right)=\frac{1}{M}{\displaystyle \sum _{j=1}^M{c}_{j,n+1}}(t)\end{array}$$

(3)

The relationship between the amplitude, times, and denoising effect of the added white Gaussian noise into the noisy signal x(t) is represented as follows:

$${\epsilon }_i=\frac{\epsilon }{\sqrt{M}}$$

(4)

where ${\epsilon }_i$ is the denoising effect of the added white Gaussian noise, and ε is the amplitude of the white Gaussian noise.

As shown in Formula (4), when the noise amplitude is fixed, the larger the number of times white Gaussian noise is added, the closer the final decomposition result is to the true values. If the amplitude of the added noise is too small, noise does not affect the selection of the extreme much, thus playing no small role in complementing the scales.

(5) Finally, the noisy signal x(t) can be decomposed as follows:

$$x\left(t\right)={\displaystyle \sum _{i=1}^n{c}_i}(t)+c_{n+1}(t)$$

(5)

2.2. Basis of Hankel Matrix and MSPCA

For a single signal, only one non-zero singular value corresponds to an SVD. Therefore, a singular value cannot be used for denoising. Hence, we constructed a p × q order Hankel matrix H_i from c_i:

$$H_i={\left(c_{i,km}\right)}_{p\times q}=\left[\begin{array}{cccc}{c}_{i,}{}_{11}& c_{i,12}& \cdots & c_{i,1q}\\ c_{i,21}& c_{i,22}& \cdots & c_{i,2q}\\ ⋮& ⋮& ⋮& ⋮\\ c_{i,p1}& c_{i,p2}& \cdots & c_{i,pq}\end{array}\right]$$

(6)

where c_i, _km = c_i (k + m − 1), and k + m − 1 = N, k and m are the row and column identities, respectively, and N is the total sampling number of c_i.

Based on the SVD, matrix H_i can be decomposed as follows: $H_i=U\sum V^{\mathrm{T}}$, ($U\in {\mathrm{R}}^{p\times p}$, $\sum \in {\mathrm{R}}^{p\times q}$ and $V\in {\mathrm{R}}^{q\times q}$). We marked $T=U\sum$ and $P=V$. H can be then treated as the sum of the outer products of the q vectors:

$$H_i={\mathrm{TP}}^{\mathrm{T}}={\displaystyle \sum _{k=1}^q{t}_k{p}_k{}^{\mathrm{T}}}=t_1{p}_1{}^{\mathrm{T}}+t_2{p}_2{}^{\mathrm{T}}+\cdots +t_q{p}_q{}^{\mathrm{T}}$$

(7)

where t_k and p_k indicate the score vector and load vector, respectively, k = 1, 2,…, q, and T represents the matrix transposition.

The t_k length reflects the $H_i$ coverage in direction p_k: a longer t_k represents greater importance, while a shorter t_k is mainly noise and placed in E. Subsequently, matrix H can be rewritten as

$$H_i=t_1^{\prime }{p}_1{}^{\mathrm{T}}+t_2^{\prime }{p}_2{}^{\mathrm{T}}+\dots +t_{q^{\prime }}^{\prime }{p}_{q^{\prime }}{}^{\mathrm{T}}+E$$

(8)

where $t_k^{\prime }$(k = 1, 2, ……, q’) is the score vector; a larger t_k remains unaltered, whereas a shorter t_k is set as a zero vector.

Considering that the PCA of the matrix, $H_i$ is equivalent to the eigenvector analysis of the covariance matrix $H_i{}^{\mathrm{T}}{H}_i$, and the load vector of matrix $H_i$ is the same as the eigenvector of $H_i{}^{\mathrm{T}}{H}_i$. Therefore, we determine $t_j^{\prime }$ depending on the eigenvalues of the matrix $H_i{}^{\mathrm{T}}{H}_i$, and the descending ordered eigenvalues are

$${\widehat{\lambda }}_1\ge {\widehat{\lambda }}_2\ge \cdots \ge {\widehat{\lambda }}_q$$

(9)

In general, 85% of the cumulative ratio between an eigenvalue and the total sum of the eigenvalues is selected as the threshold to determine $t_k^{\prime }$. The eigenvalues corresponding to the abandoned eigenvectors are set as zero vectors, thus providing the reconstructed Hankel matrix as

$${\widehat{H}}_i=\widehat{\mathit{T}}{P}^{\mathrm{T}}=t_1^{\prime }{p}_1{}^{\mathrm{T}}+t_2^{\prime }{p}_2{}^{\mathrm{T}}+\dots +t_{q^{\prime }}^{\prime }{p}_{q^{\prime }}{}^{\mathrm{T}}$$

(10)

The PCA-denoised component $c_i^{\prime }$ is obtained from the data of the first row and last column in ${\widehat{H}}_i$.

2.3. EEMD-MSPCA-Based Denoising Method

The first few components obtained by EEMD are mainly high-frequency noises, and noise may exist in the low-frequency components. Therefore, we proposed an EEMD-MSPCA-based seismic signal denoising method (Figure 1), and the specific steps are as follows.

(1) The noisy signal x(t) is decomposed by EEMD described in Section 2.1, obtaining n IMFs $c_i(t)$ and one residual component $c_{n+1}(t)$.

(2) The variance contribution rates (VCR) of the IMFs are calculated, and the first few high-frequency components with VCRs < 0.01 are removed. The VCR of the ith IMF $c_i(t)$ is defined as

$$VC{R}_i=\left[\frac{1}{N}{\displaystyle \sum _{t=1}^N{c}_i{(t)}^2}-{\left(\frac{1}{N}{\displaystyle \sum _{t=1}^N{c}_i(t)}\right)}^2\right]/{\displaystyle \sum _{i=1}^{n+1}\left[\frac{1}{N}{\displaystyle \sum _{t=1}^N{c}_i{(t)}^2}-{\left(\frac{1}{N}{\displaystyle \sum _{t=1}^N{c}_i(t)}\right)}^2\right]}$$

(11)

(3) The Hankel matrix $H_i$ is constructed for each remaining IMF $c_i(t)$, which is treated as the data matrix for PCA denoising, and its score vector and load vector are calculated through SVD.

(4) Eigenvalues and eigenvectors of the covariance matrix $H_i{}^{\mathrm{T}}{H}_i$ are calculated, and the eigenvalues are ranked in descending order. The score vector that corresponds to an eigenvalue with a cumulative eigenvalue percentage larger than 85% is set as a zero vector, and a detailed illustration is shown in Section 3. Subsequently, the Hankel matrix is reconstructed with the new score vector and load vector to obtain the denoised signal $c_i^{\prime }(t)$.

(5) The signal $c_i^{\prime }(t)$ is further denoised with soft thresholding and calculated as follows:

$$T{h}_{\mathrm{s}}=\{\begin{array}{l}\mathrm{sign}(c_i^{\prime })(\left|c_i^{\prime }\right|-T),\left|c_i^{\prime }\right|>T\\ 0,\left|c_i^{\prime }\right|\le T\end{array}$$

(12)

where |·|, sign(·), and $T$ indicate the absolute value, sign function, and threshold value, respectively. A fixed threshold $T=\sigma \cdot \sqrt{2\log (N)}$ ($\sigma$ denotes data variance) was used in this study.

(6) For the remaining IMFs, denoised components ${\widehat{c}}_i(t)$ can be obtained by repeating steps (3–5), and the final denoised signal $\widehat{x}(t)$ is obtained by stacking the denoised components ${\widehat{c}}_i(t)$, as follows:

$$\widehat{x}(t)={\widehat{c}}_i(t)+\cdots +{\widehat{c}}_n(t)+{\widehat{c}}_{n+1}(t)$$

(13)

3. Synthetic Tests of EEMD-MSPCA-Based Denoising

Blocks, ECG, Bumps, Doppler shift, and Heavy sine signals with white Gaussian or random noise are usually adopted to test a denoising method, because it is easy to calculate the quantitative parameter signal to noise ratio (SNR) compared with a real signal. Thus, it is convenient for other researchers to reproduce or compare their results with ours. The Blocks signal [s(t) in Figure 2a] is selected to demonstrate the proposed denoising algorithm. The number of signal sampling points was set as N = 1024, and the added noise was expressed as n(t) = 0.2 × randn(1, 1024), i.e., the noise is 0.2 times the standard normal distribution, and the noisy signal x(t) is shown in Figure 2b.

With EEMD, x(t) is decomposed into nine IMF components, c₁–c₉, and one residual component, c₁₀ (Figure 3a). These data indicate that the first few IMF components are mainly high-frequency noises, and the results obtained via Formula (11) show that the VCRs of the first two IMFs are <0.01. Therefore, they are eliminated (marked as zero vectors in Figure 3b) before PCA denoising.

c₃ was selected to illustrate the PCA-based denoising process. First, the Hankel matrix $H_3$ constructed with c₃ using Formula (6) is decomposed into $H_3=U\sum V^{\mathrm{T}}$ through SVD, with $T=U\sum$ and $P=V$. Eigenvalues and eigenvectors of the covariance matrix $H_3{}^{\mathrm{T}}{H}_3$ are calculated, and the eigenvalues are ranked in descending order (red points in Figure 4). The percentages of cumulative eigenvalues are shown in blue points in Figure 4. The score vectors corresponding to the cumulative eigenvalue percentages <85% are used for PCA-based c₃ denoising, and the denoised signal ${c^{\prime }}_3$ is obtained by the reconstructed Hankel matrix ${\widehat{H}}_3$, using Formula (10). Finally, the soft thresholding-denoised ${\widehat{c}}_3$, derived from Formula (12), is shown in Figure 3b. By utilising the same process, the remaining components are denoised to obtain ${\widehat{c}}_4~{\widehat{c}}_{10}$.

As illustrated in Figure 3b, the denoised ${\widehat{c}}_3$ enhanced the local amplitude characteristics of N = 100–250 and N = 300–450 and suppressed the noise amplitudes of other bands to an extent. Similar denoising results were obtained for the other IMFs, demonstrating the effectiveness of the EEMD-MSPCA-based denoising method. The final denoised signal (Figure 2c) was obtained using $\widehat{x}(t)={\displaystyle \sum _{i=3}^9{\widehat{c}}_i(t)+{\widehat{c}}_{10}(t)}$. SNR was selected to quantitatively evaluate the denoising performance, as defined in Formula (14). The SNRs increased from 6.99 to 12.55, and the denoised signal effectively retained the signal characteristics.

$$SNR=10\lg [{\displaystyle \sum _{i=1}^N{s}_i^2/}{\displaystyle \sum _{i=1}^N{n}_i^2}]$$

(14)

where $s_i$ is the amplitude of the signal without noise and $n_i$ is the amplitude of noise, which is equal to the difference between $s_i$ and the denoised signal.

Furthermore, the ECG, Bumps, Doppler shift, and Heavy sine signals (separately shown top in Figure 5a–d) were selected to test the proposed denoising method. The noise is added using the same steps as those used in the Blocks signal test. The obtained noisy signals and EEMD-MSPCA-denoised signals are shown in the middle and bottom of Figure 5a–d, respectively, and their SNRs are listed in

Table 1. SNRs of the five noisy signals with different denoising methods.

Signal Name	Noisy Signal	Low-Pass Filter	Wavelet Reconstruction	EEMD Reconstruction	Hankel–SVD	EEMD-Hankel–SVD	Wavelet-MSPCA	EEMD-MSPCA
Blocks	6.99	4.67	10.32	10.45	8.69	10.78	11.83	12.55
ECG	0.49	4.46	8.32	8.63	2.93	9.68	9.26	9.43
Bumps	12.55	6.29	17.83	18.38	14.93	19.23	19.88	20.13
Doppler shift	9.32	6.16	15.18	15.59	11.46	16.51	16.58	16.84
Heavy sine	9.51	15.95	18.52	18.67	14.87	19.18	19.04	19.28

. Figure 5 depicts the local characteristics of the denoised signals, which are obviously enhanced, and Table 1 displays the increased SNRs of the denoised signals. A comparison of the proposed method with other denoising methods will be shown in the Discussion section.

4. MS Signal EEMD-MSPCA Denoising and P-Phase Arrival Picking

To further verify the effectiveness of the EEMD-MSPCA-based denoising method, 500 MS signals with a sampling frequency of 6000 Hz were randomly selected from the MS monitoring system in the Yongshaba mine (China) as the application dataset, which usually contains different noises (e.g., rock drilling, loco transportation, electrical noise, mine drawing noise, and blower vibration noise) and acts as a good dataset to test the denoising method. Because the real MS signal is unknown, we used the automatic P-phase arrival picking error instead of the signal SNR to evaluate the denoising performance. The STA/LTA, PAI-K, and AIC picking methods described previously [25] were employed, and the picking error was defined as the difference between the automatic and manual picking points; a smaller picking error represents better denoising performance.

Three typical MS signals before and after denoising and the P-phase arrival picking results are shown in Figure 6. The original waveforms (Figure 6a–c) demonstrate that there is heavy noise in the recorded signals with low SNRs, especially for the third signal (Figure 6c), which includes high-frequency noise, power frequency noise, and spike noise. The EEMD-MSPCA-based denoising method can effectively remove the above-mentioned noises and retain the characteristics of P-phase arrival, which is favourable for P-phase arrival picking. The STA/LTA, kurtosis, and AIC values of the denoised signals are sharper than those without denoising, verifying the improvement of signal SNRs and the effectiveness of the proposed denoising method.

The picking error statistics of the STA/LTA picking, PAI-K picking and AIC picking for the 500 MS signals are shown in Figure 7a–c, respectively. The data indicate that there are more small picking error signals after the EEMD-MSPCA denoising, and the signal number of picking errors that are larger than 30 ms or no pickings is reduced, suggesting an improvement in SNR by the proposed denoising method. Furthermore, the total cost function (TCF) proposed by Li et al. [26] was selected for quantitative comparison. The cost function (CF_i) and TCF are defined in Formulas (15) and (16), respectively, where the TCF is obtained through the summation of cost functions corresponding to P-phase arrival picking errors, and a small cost is incurred for a small picking error. Therefore, a smaller TCF indicates better picking results. Figure 7 shows that the TCFs of the denoised signals are smaller than those based on the recorded signals, and the TCF based on denoised signals and the AIC method showed the best picking results.

$$C{F}_i=\{\begin{array}{cc}0.0,& \mathrm{accuracy}<5\mathrm{ms}\\ 0.2,& 5\mathrm{ms}\le \mathrm{accuracy}<10\mathrm{ms}\\ 0.4,& 10\mathrm{ms}\le \mathrm{accuracy}<15\mathrm{ms}\\ 0.6,& 15\mathrm{ms}\le \mathrm{accuracy}<20\mathrm{ms}\\ 0.8,& 20\mathrm{ms}\le \mathrm{accuracy}<25\mathrm{ms}\\ 1.0,& 25\mathrm{ms}\le \mathrm{accuracy}<30\mathrm{ms}\\ 1.5,& \mathrm{accuracy}\ge 30\mathrm{ms}\end{array}$$

(15)

$$TCF={\displaystyle \sum _{i}C{F}_i}$$

(16)

5. Discussion

To better validate the proposed denoising method, low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, EEMD-Hankel–SVD, and wavelet-MSPCA-based denoising methods were selected for comparison. The filter used is a four-order low-pass Butterworth filter with cutoff frequency of 200 Hz, wavelet reconstruction and EEMD reconstruction both remove the first few detail coefficients with VCRs smaller than 0.01, Hankel–SVD was proposed by Jiang et al. [23], EEMD-Hankel–SVD was proposed by Li et al. [27], and the wavelet-MSPCA-based denoising was proposed by Bakshi [13]. Synthetic tests were performed on a Thinkpad laptop with 8 GB RAM and a 1.6 GHz Intel Core i5 Dual Core Processor. The results showed that the computational costs of the low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, EEMD-Hankel–SVD, wavelet-MSPCA, and EEMD-MSPCA-based denoising methods were 0.01, 0.02, 0.25, 0.18, 0.15, 2.96, and 3.03 s, respectively. The EEMD has a relatively large computational cost because of the addition of white Gaussian noise in each EMD step. The EEMD-Hankel–SVD and proposed EEMD-MSPCA-based denoising methods require a much longer computation time than wavelet-MSPCA-based denoising, because the EEMD/wavelet decomposition times are proportional to the row number of the Hankel matrix shown in Formula (6). The denoising performance of the seven methods for the noisy Blocks signal (Figure 2b) is shown in Figure 8, and the SNRs of different denoising methods for the five noisy signals used in the synthetic tests are listed in Table 1.

Figure 8 demonstrates that the low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, and wavelet-MSPCA-based denoising methods are inferior to the EEMD-MSPCA-based denoising method in some respects. Moreover, Table 1 shows that The SNR of the ECG signal improved from 0.49 to >8. Figure 5a shows that many amplitudes in the time series of the ECG signal are zero, causing the ${\displaystyle \sum _{i=1}^N{s}_i^2}$ and ${\displaystyle \sum _{i=1}^N{n}_i^2}$ in the SNR definition to be relatively small and large, respectively, after including noises, and a small SNR can then be obtained. Furthermore, there is a coefficient of 10 for the SNR definition. Thus, denoising can significantly improve the SNR in this case. Thus, the low-pass filter improves the waveform quality; however, it can decrease the SNR due to waveform shift, which is harmful for P-wave arrival picking. All the other denoising methods improved the SNRs, with the EEMD-Hankel–SVD and EEMD-MSPCA obtaining the best denoising performance. The underlying reasons are as follows: Hankel–SVD denoising removes information corresponding to low eigenvalues, while the other denoising methods eliminate high-frequency components in the denoising process, which illustrates the importance of removing high-frequency noise. Subsequently, EEMD reconstruction can reduce the mode mixing effect with a slightly superior denoising result compared to the wavelet reconstruction method. Finally, PCA- and soft thresholding-based denoising can denoise the remaining low-frequency noise [28], as well as retain the signal characteristics.

6. Conclusions

In this study, an EEMD-MSPCA-based signal denoising method was proposed to overcome the difficulty of wavelet-MSPCA in denoising a single signal. This method was verified using synthetic tests and MS signal denoising. The main conclusions are summarised as follows: (1) The proposed method eliminates high-frequency noise through EEMD and the variance contribution rate; subsequently, a Hankel matrix is constructed with each IMF for PCA denoising. PCA and soft thresholding can denoise low-frequency components, as well as retain signal characteristics. (2) EEMD-MSPCA-based denoising performs better than the low-pass filter, wavelet reconstruction, EEMD reconstruction, Hankel–SVD, EEMD-Hankel–SVD, and wavelet-MSPCA-based denoising methods. (3) EEMD-MSPCA-denoised signals effectively improve the P-phase arrival picking accuracy of the STA/LTA, PAI-K, and AIC methods, and the EEMD-MSPCA and AIC picking method achieves a good P-phase arrival picking result, which plays an important role in the MS parameter calculation (e.g., source location and focal mechanism). Further denoising improvements can consider other adaptive decomposition methods (e.g., VMD, local mean decomposition (LMD)), noisy component selection parameters (e.g., signal entropy and coefficient), and combining with other denoising methods.