Research on De-Noising Method of Grounded Electrical Source Airborne Transient Electromagnetic Data Based on Singular Spectrum Analysis

Abstract

The grounded electrical source airborne transient electromagnetic (GREATEM) system is widely used in groundwater resources detection, geothermal resource detection, geological structure detection, and other fields due to its wide detection range, high detection efficiency, and high resolution. The field data received by the GREATEM system is easily affected by various noises, such as instrument system noise, power frequency noise, sferics noise, and other noise, which reduce the data signal-to-noise ratio (SNR) and affects the data interpretation accuracy. This paper proposes a singular spectrum analysis (SSA) for the GREATEM data de-noising in response to this problem. First, we calculate the electromagnetic response of a uniform half-space using a GREATEM system with an electrical source to verify the effectiveness of the SSA algorithm for GREATEM data de-noising. To determine the appropriate parameters for SSA, we propose a particle swarm optimization algorithm to choose the window length. Later, SSA is used to decompose a synthetic quasi-two-dimensional earth model of GREATEM data. After SSA, the SNR of the reconstructed signal increased by 36 dB, and the RMSE does not exceed 4.9 × 10⁻⁶, which verifies the feasibility of the SSA for de-noising GREATEM data. Finally, through field measurement data processing, the effectiveness of the method is further confirmed.

1. Introduction

The grounded electrical source airborne transient electromagnetic (GREATEM) system is widely used in geophysical exploration due to its wide detection range, high detection efficiency, and high resolution [1,2,3,4,5,6,7,8]. During in- flight measurements, the electromagnetic signals received by GREATEM are easily affected by various noises, such as motion noise, electromagnetic noise, etc. Due to these noises, the signal-to-noise ratio (SNR) of the observing secondary field with weaker energy is reduced, which affects the accuracy of the data interpretation results.

There are currently many methods to eliminate a specific type of electromagnetic noise for the GREATEM data. Wang et al. [9] proposed a wavelet-based method to remove the baseline drift of GREATEM signals. Li et al. [10] proposed a comprehensive wavelet de-noising method based on a sym8 wavelet base, effectively suppressing white noise and baseline drift low-frequency noise in GREATEM signals. Li et al. [11] used stationary wavelet transformation to remove ground spike noise and background noise in the GREATEM data. Liu et al. [12] proposed the ensemble empirical mode decomposition (EEMD) method to suppress the noise caused by motion in the GREATEM system. Ji et al. [13] proposed an exponential fitting-adaptive Kalman filter (EF-AKF) based on GREATEM signals’ characteristics and achieved significant de-noising effects. Wu et al. [4] combined a wavelet and neural network to predict and filter high-frequency motion noise in ATEM signals and achieved good results. Li. et al. [14] adopted the EEMD-based adaptive filtering method, which suppresses the baseline drift in the GREATEM signal. Wu et al. [15] proposed a deep learning method based on a de-noising autoencoder (DAE). By constructing a training set of a large amount of simulated data, this can simultaneously deal with multiple types of noise in AEM data and reduce the influence of human subjective factors. Wu et al. [16] proposed a variational Bayesian-based adaptive Kalman filter (VBAKF) to de-noise the GATEM data. The method is based on the results of an exponential fit, which is suitable for exponentially decaying signals. With the development of computer science, many researchers began to use artificial intelligence algorithms to invert electromagnetic data. He et al. [8] developed a three-dimensional inversion of semi-airborne transient electromagnetic data based on a particle swarm optimization-gradient (PSO) descent algorithm. Wu et al. [17] proposed a deep learning method to estimate the earth resistivity model for ATEM observation. These methods perform the de-noising and inversion functions at the same time.

Adaptive filtering methods require the collection of noisy data near the measured data for adaptive cancellation in the above methods. The correlation between the noise data and the original data is often too significant, resulting in poor noise elimination. The wavelet de-noising method needs to select a suitable wavelet base and threshold according to the characteristics of the noise, which is not universal [4,9,10,11]. The EMD may miss helpful information and has a modal aliasing problem due to its frequency band segmentation characteristics [12,14,18,19].

Singular spectrum analysis (SSA) was first proposed by Broomhead and King [20] and developed from the Karhumen–Lèove decomposition theory [21]. It is a new time series analysis method, which is used chiefly in signal de-noising and prediction. It shows excellent gravity and magnetic potential field separation and magnetic resonance sounding (MRS) signal de-noising [22,23]. Therefore, we proposed to use singular spectrum analysis to remove noises in GREATEM data.

2. Singular Spectrum Analysis

SSA is a method to study the de-noising or prediction of one-dimensional time-series signals. It mainly includes two essential parts: decomposition and reconstruction. Decomposition includes embedding and SVD. Reconstruction includes regrouping and diagonal averaging. The specific steps are as follows.

(1)

Embedding

One-dimensional TEM data $\mathit{y}=\left[y\left(1\right),y\left(2\right),\dots ,y\left(N\right)\right]$ can be mapped into an M × L trajectory matrix Y

$$Y=\left[\begin{array}{ccc}y\left(1\right) y\left(2\right)& \cdots \cdots & y\left(L\right)\\ y\left(2\right) y\left(3\right)& \cdots \cdots & y\left(L+1\right)\\ ⋮ ⋮& \cdots \cdots & ⋮\\ y\left(M\right) y\left(M+1\right)& \cdots \cdots & y\left(N\right)\end{array}\right]$$

(1)

where M is the window length and L = N − M + 1, the trajectory matrix Y is called the delay matrix [24].

(2)

Singular Value Decomposition

First, calculate the covariance matrix C = YY^T of the trajectory matrix Y, which is expressed as

$$C=\left[\begin{array}{ccc}c\left(0\right) c\left(1\right)& \dots & c\left(M-1\right)\\ c\left(1\right) c\left(0\right)& \dots & c\left(M-2\right)\\ ⋮ ⋮& \ddots & ⋮\\ c\left(M-1\right) c\left(M-2\right)& \dots & c\left(0\right)\end{array}\right]$$

(2)

C is a Toeplitz matrix and needs to obtain its eigenvalue λ and eigenvector matrix U. The eigenvalues are arranged in descending order and correspond to M × 1 eigenvector. The eigenvalue vector is expressed as

$$\mathsf{\lambda }=\left[\begin{array}{c}{\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2,\cdots ,{\mathsf{\lambda }}_M\end{array}\right], \left|\begin{array}{c}\begin{array}{c}{\mathsf{\lambda }}_1\end{array}\end{array}\right|\ge \left|\begin{array}{c}{\mathsf{\lambda }}_2\end{array}\right|\ge \cdots \ge \left|\begin{array}{c}{\mathsf{\lambda }}_M\end{array}\right|$$

(3)

The singular value vector σ of the matrix C can be presented as

$$\mathsf{\sigma }=\left[\begin{array}{c}\begin{array}{c}{\mathsf{\sigma }}_1,{\mathsf{\sigma }}_2,\cdots ,{\mathsf{\sigma }}_M\end{array}\end{array}\right]=\left(\sqrt{\left|\begin{array}{c}{\mathsf{\lambda }}_1\end{array}\right|}, \sqrt{\left|{\mathsf{\lambda }}_2\right|}\cdots ,\sqrt{\left|\begin{array}{c}{\mathsf{\lambda }}_M\end{array}\right|}\right)$$

(4)

In general, the singular value vector σ of matrix C is called the singular spectrum of Y. The effective singular value is the singular value that is greater than zero in the singular value vector. In the actual data processing, the elements in the singular spectrum vector σ of the real sampled signal vector are generally larger than zero. The part with the smaller singular value is close to zero. The part that is more significant than zero is usually regarded as an effective singular value.

The trajectory matrix Y is composed of several elementary matrices:

Y = Y₁ + Y₂ + …+Y_J

(5)

Let $=U_j/\sqrt{{\lambda }_j}$, where $U_j$ and $V_j$ are the left and right eigenvectors corresponding to the eigenvalues after the trajectory matrix Y is decomposed. Then the trajectory matrix of the jth component in Y can be expressed as

$$Y_j=\sqrt{{\lambda }_j}{U}_j{V}_j^T j=1,2,\cdots ,M$$

(6)

Substituting ${\mathbf{V}}_j$ in (6) and merging similar terms, ${\mathbf{Y}}_j$ can be presented as

$$Y_j=U_j{U}_j^{T}Y j=1,2,\cdots ,M$$

(7)

(3)

Reconstruction

Let $J=\left\{j_1,j_2,\cdots ,j_d\right\}$ be the indices corresponding to d eigenvalues, then the matrix $Y_J$ will be synthesized as

$$Y_{\mathrm{J}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{d}}{Y}_{{\mathrm{j}}_{\mathrm{i}}}$$

(8)

Divide the subscripts of the matrix $Y_j$ into p disjoint subsets $J_1,J_2,\cdots J_p$, then the original trajectory matrix Y can be expressed as

$$Y={{\displaystyle \sum }}_{j=1}^p{Y}_{J_p}$$

(9)

(4)

Diagonal averaging

Finally, the trajectory matrix is converted into one-dimensional data with N time channels in the diagonal averaging operation. $\widehat{S}$ is the L × K trajectory matrix. According to the diagonal averaging formula, $\widehat{S}\left(n\right)$ is converted into a one-dimensional time series. The diagonal averaging formula is as

$$\widehat{S}\left(n\right)=\{\begin{array}{c}\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\widehat{s}}_{i,n-i+1} for 1\le n<M\\ \frac{1}{M}{\displaystyle {\displaystyle \sum }_{i=1}^M}{\widehat{s}}_{i,n-i+1} for M\le n<L\\ \frac{1}{N-n+1}{\displaystyle {\displaystyle \sum }_{i=n-L+1}^{N-L+1}}{\widehat{s}}_{i,n-i+1} for L\le n<N\end{array}$$

(10)

Different estimated signal $s_i$ corresponds to the corresponding singular values ${\mathsf{\sigma }}_i$, where the larger the singular value ${\mathsf{\sigma }}_i$ is, the closer the signal is to the expected signal after de-noising.

3. Simulation Analysis Using SSA for GREATEM Signal De-Noising

In order to verify the effectiveness and correctness of the SSA algorithm for GREATEM detection data de-noising, we calculated the electromagnetic response of a uniform half-space using a GREATEM system with an electrical source. The specific calculation parameters are that the length of the long wire is L = 1000 m, the transmission current is I = 10 A, the effective area of the receiving coil is S = 10,000 m², and the time series length of the synthetic signal is 5001. For the simulation calculation process, refer to reference Ji et al. [13]. Simulation of GREATEM signals: a long grounded wire source is used, the coordinate origin is set in the middle of the wire, the length of the half wire is L, and the vertical component of the induced electromotive force is expressed in the layered earth model as Equation (11) shown

$$V=-i\omega I\frac{{\mu }_{0}S}{4\pi }{{\displaystyle \int }}_{-L}^L\frac{y}{R_d}{{\displaystyle \int }}_0^{\infty }\left(1+r_{TE}\right)e^{u_{0}z}\frac{{\lambda }^2}{u_0}{J}_1\left(\lambda R\right)d\lambda d{x}^{\prime }$$

(11)

where i is an imaginary unit, where f is the frequency (Hz), I is the emission current, representing the air permeability, S is the equivalent area of the receiving coil, $R_d={\left[{\left(x-x^{\prime }\right)}^2+y^2\right]}^{\frac{1}{2}}$ is the distance between the transmitting source and the receiving coil, and λ is the Hankel transform the integral variable. r_TE is the reflection coefficient, and J₁ is the first-order Bessel function. The resistivity of uniform half-space is 100 Ω·m.

Calculate the electromagnetic response of the uniform half-space of the electrical source, and set the simulation parameters as follows: the emission current I = 10 A, the source length of the half-length wire is 500 m, the equivalent area of the receiving coil is S = 10,000 m², and the ideal ground-air transient electromagnetic signal curve can be obtained, whose time track length is 5001. Figure 1 shows the calculated GREATEM response of uniform half-space.

3.1. Parameter Selection

The singular spectrum analysis method has two parameters: the window length M and the effective singular value order. Whether the user can appropriately select two parameters is the key factor to determining the de-noising effect of the SSA for GREATEM data.

3.1.1. Window Length

The first step in the SSA algorithm is to select the appropriate window length M to embed in the time trajectory matrix Y. The value of the window length M directly determines the de-noising effect. Generally speaking, the window length M should be less than half of the trajectory matrix length, which is determined by the properties of the matrix. Figure 2 shows the comparison result of the de-noising effect on the calculated theoretical GREATEM response of uniform half-space with power frequency noise, where the window length M is the value of 30, 350, 1500, and 2500, respectively.

As shown in Figure 2, for the length of 5001 synthetic signals, the window length 30 is too small, the noise reduction effect is not ideal, and the signal still contains apparent noise. The window length 1500 or 2500 is too large, the separated signal decay ends prematurely, and the useful signal components are not properly retained.

In applying the SSA method, some scholars usually take the window as half of the length of the trajectory matrix Y [25,26] and have achieved good results. Figure 2 concludes that this window selection method could not appropriately separate GREATEM signals and noise. Therefore, we study the window length selection in the de-noising of GREATEM data using the SSA algorithm. We select the window length M value interval in [10, 1000], the step length is 10, and use root means square error (RMSE) as the evaluation index. The smaller the RMSE, the better the separation effect. Figure 3 shows the signal separation results for different window lengths. We can conclude that the appropriate window length is not a specific value but an interval. Select the window length in this interval, and the separated GREATEM signal effect is similar.

Therefore, we propose to use a particle swarm optimization (PSO) algorithm to select the appropriate window length M more accurately. The PSO algorithm is a kind of self-adaptive random optimization technology which has the characteristics of a simple iterative format and fast convergence to the optimal solution of the equation. The basic principle of the PSO algorithm is:

(1)

Let m particles form a group in an S-dimensional search space, and the ith particle is represented as an S-dimensional vector ${\mathit{X}}_i=\left(x_{i1},x_{i2},\dots ,x_{iS}\right),i=1,2,\dots ,m$; then we substitute ${\mathit{X}}_i$ into the objective function to calculate the corresponding fitness value.

(2)

Let the optimal individual extremum of the current iteration of the ith particle be ${\mathit{P}}_i=\left(p_{i1},p_{i2},\dots ,p_{iS}\right)$; the velocity is ${\mathit{V}}_i=\left(v_{i1},v_{i2},\dots ,v_{iS}\right)$; the optimal global extremum of the particle swarm current iteration is ${\mathit{P}}_g=\left(p_{g1},p_{g2},\dots ,p_{gS}\right)$.

(3)

Update the speed and position of the ith particle in S dimension according to Equations (12) and (13).

$${\mathit{V}}_i\left(t+1\right)=\omega {\mathit{V}}_i\left(t\right)+c_1{v}_1\left({\mathit{P}}_i\left(t\right)-{\mathit{X}}_i\left(t\right)\right)+c_2{v}_2\left({\mathit{P}}_g\left(t\right)-{\mathit{X}}_i\left(t\right)\right)$$

(12)

$${\mathit{X}}_i\left(t+1\right)={\mathit{X}}_i\left(t\right)+{\mathit{V}}_i\left(t+1\right)$$

(13)

In Equations (12) and (13), t is the number of iterations; ω is the inertia weight, generally take $\omega ϵ \left[0.4,1.0\right]$; $c_1$ and $c_2$ are learning factors, generally take $c_1,c_2 ϵ \left[0, 2\right]$;$r_1$ and $r_2$ are mutually independent and obey a pseudo-random number of the uniform distribution on [0, 1].

Set the parameters to ω = 0.7, $c_1=c_2=2$, and the population size to 20. At the same time, the initial population position range is limited [10, 1000] [27,28]. Using permutation entropy as a fitness function of particle swarm optimization, the window length and fitness function curve shown in Figure 4 is obtained through 20 PSO algorithm iterations. As shown in Figure 4, the algorithm converges at 155, which is similar to Figure 3. It proves the effectiveness of particle swarm optimization.

3.1.2. Reconstruction Order

Many scholars use dichotomy, the mean truncation method, etc., to reconstruct the signal in selecting the reconstruction order of the SSA algorithm [26]. To choose the reconstruction order for GREATEM data de-noising, we use noisy simulation and measured data to study. Figure 5 shows three groups of different GREATEM signals and their singular spectra. Figure 5a is a calculated GREATEM signal with white noise; Figure 5b is a calculated GREATEM signal with power frequency noise and a little white noise; Figure 5c is a set of GREATEM attenuation curves measured in front of the Geological Palace in Changchun. We can conclude that the singular values of GREATEM signals are the first two. Therefore, the SSA algorithm is used for GREATEM data de-noising, and the reconstruction order is 2.

3.2. De-Noising Simulation Analysis for GREATEM Signal with Measured Noise

In the process of GREATEM system detection, the measured data often contains different types of noise. When the detection is near a transformer or under high-voltage transmission lines, the noise is mainly caused by the higher harmonic noise. The attenuation law of GREATEM signals is close to the exponent function. The early attenuation rate of the signal is fast, and its amplitude is large. The late measured data is similar to DC data, and its amplitude is 4–5 orders lower than that of the early period. Since the late data is weak, the noise interference of the late-stage signal is more serious. Firstly, we simulated and calculate the GREATEM signal with measured noise, as shown in Figure 6. We can see that the effect of measured noise on the late signal is more significant. Then we used SSA to de-noise the signal, as shown in Figure 7. The GREATEM signal is reconstructed according to the different singular values, and the noise components are separated from it.

Figure 7a is the reconstructed GREATEM signal, and Figure 7b is the separated measured noise. Figure 8 shows the signal comparison before and after SSA de-noising. The Xn is the calculated GREATEM signal with noise, the SSA is the reconstructed GREATEM signal by SSA, and Sn is the calculated GREATEM signal without noise. After comparison with Figure 6, the reconstructed signal fits well with the calculated GREATEM signal. This verifies that the SSA algorithm could effectively separate the GREATEM simulation signal.

From the change of the signal curve, the reconstructed signal is smooth and fits well with the original signal. We quantify the de-noising effect using SNR and RMSE. The expressions are as follows:

$$SNR=10\times {\log }_{10}\frac{P_{signal}}{P_{noise}}$$

(14)

$$RMSE=\sqrt{({{\displaystyle \sum }}_{n-1}^N{\left[\widehat{V_z}\left(n\right)-\widehat{V_z}\left(n\right)\right]}^2])/N}$$

(15)

According to the formula (14) and (15), the SNR of the noisy synthetic signal is 34.8265 dB. After SSA, the SNR is increased to 65.8314 dB, adding nearly 31 dB. The RMSE before and after SSA are 3.6335 × 10⁻⁴ and 1.0002 × 10⁻⁵, respectively, verifying that the SSA algorithm has a significant de-noising effect on GREATEM signals.

3.3. De-Noising Simulation Analysis for Quasi-Two-Dimensional Earth Model

In order to further explore the de-noising effect of the SSA algorithm, we used a quasi-two-dimensional earth model for research. The theoretical model is shown in Figure 9a, the surrounding rock resistivity was 100 Ω·m; the resistivity of the low-resistance layer was 10 Ω·m, the depth was from 150 m to 220 m, and the simulation data inversion result is shown in Figure 9b. After adding the higher harmonic noise and a little white noise to the simulation data, before and after de-noising were inverted. Figure 9c is the inversion result before de-noising. Due to the higher harmonic noise and the white noise, the inversion result can not reflect the actual model, and the position and thickness of the low-resistance layer have large deviations after inversion. Figure 9d is the inversion result after de-noising by the SSA algorithm, consistent with the simulation data inversion result in Figure 9d, proving the method’s effectiveness.

3.4. Comparison of SSA and Other De-Noising Methods

Many scholars use the wavelet threshold method (WT), empirical mode decomposition (EMD), and variational mode decomposition (VMD) to de-noise the GREATEM signal. We compare these four algorithms’ de-noising results to verify the superiority of SSA in de-noising GREATEM signals. We use calculated GREATEM signals of the quasi-two-dimensional earth model in Figure 9. We adopt the sym5 wavelet basis, five-layer decomposition, and the principle of fixed threshold (sqtwolog) for WT. The analysis results are shown in Figure 10.

Figure 10 shows the calculated GREATEM signal with the higher harmonic noise and the attenuation curve obtained after different de-noising methods. The Xn curve represents the synthetic GREATEM signal with noise; the Sn curve represents the calculated GREATEM signal; the WT, EMD, VMD, and SSA curves represent the de-noising results using WT, EMD, VMD, and SSA, respectively. From Figure 10, we know that all these four de-noising methods can effectively eliminate noise. Among them, the GREATEM signal separated by SSA is closer to the calculated GREATEM signal, proving SSA’s superiority in de-noising GREATEM signals.

Among these four methods, the WT method needs to choose the wavelet base and threshold according to the signal and noise characteristics in the actual de-noising process. The EMD algorithm needs to determine the IMF component structure rules’ weight and has end effects. SSA has fewer parameters, easy calculation, and high efficiency.

Table 1. Comparison of indexes for different algorithms to eliminate the higher harmonic noise.

Signal	SNR (dB)	RMSE
Synthetic GREATEM signal	26.0822	0.0011
Signal de-noised by WT	39.6726	2.2990 × 10−4
Signal de-noised by EMD	26.9985	9.8755 × 10−4
Signal de-noised by VMD	30.1032	6.8469 × 10−4
Signal de-noised by SSA	54.2960	4.2744 × 10−5

shows the comparison of the SNR and RMSE of different methods.

To further analyze the de-noising effect of the SSA algorithm on the measured GREATEM data, we added the measured noise to the synthetic GREATEM signal and compared four de-noising results in Figure 11. The measured noise contains white noise, spike noise, and power frequency noise. From Figure 11, we see that the signal after de-noising by the SSA is the closest to the original signal, proving that the SSA algorithm can de-noise the measured GREATEM data.

Table 2. Comparison of indexes of different methods to eliminate measured noise.

Signal	SNR (dB)	RMSE
Synthetic GREATEM signal	29.3049	0.0031
Signal de-noised by WT	51.8377	2.3079 × 10−4
Signal de-noised by EMD	49.3527	3.0723 × 10−4
Signal de-noised by VMD	47.5233	3.7892 × 10−4
Signal de-noised by SSA	64.5133	5.3631 × 10−5

shows the comparison of the SNR and RMSE of different methods. The effect of the SSA algorithm is more evident when the virtual environment contains considerable higher harmonic noise.

4. Performance of SSA with GREATEM Field Data

In November 2015, we used the GREATEM system for field detection experiments in the Daxing’an Mountain range, Heilongjiang Province. The system’s working parameters were as follows: the transmitting current was 50 A, the receiving coil area was 1080 m², the half-wire length was 744 m. The field exploration experiment site is in the base rock uplift area of the Yishu Basin. The low-resistance area is composed of sand and gravel, and the high-resistance area is horn shale. The lower resistivity part is the Yitong Basin, and the higher resistivity part is the bedrock uplift area.

We took two measured points as an example. After superimposing the measured attenuation curves, we used WT, EMD, VMD, and SSA to de-noise the attenuation curve. The results are shown in Figure 12.

Through comparison, we can see that using SSA can more effectively eliminate considerable power frequency noise, which verifies the effectiveness of the SSA method in eliminating noise in GREATEM field data and improves signal performance.

Figure 13 shows the resistivity inversion results at a depth of 150 m in Fei Hu Mountain Area, Daxing’an Mountain range measured data. Figure 13a is the investigation result of the survey area provided by Heilongjiang Institute of Geophysical and Geochemical Exploration. As shown in Figure 13b, because the signal is submerged by noise in the late stage, the low resistivity region D1 is not apparent, not effectively reflecting the real underground geological results. From Figure 13c, we can see that using SSA can effectively eliminate the false anomalies, clearly demonstrating the distribution of low-resistivity regions, and improving the imaging results’ resolution. This further verifies the effectiveness of using SSA de-noising GREATEM field data.

5. Conclusions

As a principal component analysis method of time series, SSA is mainly used in signal de-noising and prediction. When the GREATEM system is carried out on high-voltage lines or near urban areas, the data contains complex noise, which seriously affects the accuracy of the inversion interpretation results. To solve the problem, the SSA method is proposed to eliminate noise. The selection of the two critical parameters of window length and reconstruction order that affect the algorithm’s performance was analyzed based on the characteristics of GREATEM data, and using PSO, the adaptive selection of parameters was realized. Through the simulation analysis of a quasi-two-dimensional earth model, the effectiveness of the SSA method in eliminating GREATEM signal noise was verified. Through the comparison with the classic de-noising method, the superiority of SSA in de-noising GREATEM data was demonstrated, and the SSA algorithm was used to process the measured data in the Daxing’an Mountain range area to verify the practicability of the method.