Open AccessArticle

De-Noising of Magnetotelluric Signals by Discrete Wavelet Transform and SVD Decomposition

College of Geoexploration Science and Technology, Jilin University, Changchun 130026, China

Key Laboratory of Applied Geophysics, Ministry of Natural Resources, Changchun 130026, China

Changbai Mountain Volcano Comprehensive Geophysics Field Science Observation and Research Station of Ministry of Education, Changchun 130026, China

Author to whom correspondence should be addressed.

Remote Sens. 2021, 13(23), 4932; https://doi.org/10.3390/rs13234932

Submission received: 27 September 2021 / Revised: 22 November 2021 / Accepted: 1 December 2021 / Published: 4 December 2021

(This article belongs to the Topic High-Resolution Earth Observation Systems, Technologies, and Applications)

Download

Browse Figures

Review Reports Versions Notes

Abstract

Magnetotelluric (MT) sounding data can easily be damaged by various types of noise, especially in industrial areas, where the quality of measured data is poor. Most traditional de-noising methods are ineffective to the low signal-to-noise ratio of data. To solve the above problem, we propose the use of a de-noising method for the detection of noise in MT data based on discrete wavelet transform and singular value decomposition (SVD), with multiscale dispersion entropy and phase space reconstruction carried out for pretreatment. No “over processing” takes place in the proposed method. Compared with wavelet transform and SVD decomposition in synthetic tests, the proposed method removes the profile of noise more completely, including large-scale noise and impulse noise. For high levels or low levels of noise, the proposed method can increase the signal-to-noise ratio of data more obviously. Moreover, application to the field MT data can prove the performance of the proposed method. The proposed method is a feasible method for the elimination of various noise types and can improve MT data with high noise levels, obtaining a recovery in the response. It can improve abrupt points and distortion in MT response curves more effectively than the robust method can.

Keywords:

magnetotelluric; multiscale dispersion entropy; phase space reconstruction; wavelet transform; SVD decomposition

Graphical Abstract

1. Introduction

Magnetotelluric (MT) sounding can be used to explore subsurface structures by observing the natural orthogonal electric and magnetic fields. Compared with controlled-source signals, natural MT signals have the advantages of a wide frequency range, deep exploration depth, and easy-to-use equipment. Therefore, MT is widely used in many fields, such as marine exploration, volcano monitoring, and geodynamic interpretation [1,2,3]. However, during the period of observation, MT data are vulnerable to distortions from electromagnetic wave emitted by high-voltage lines, communication radio stations, and underground mining machines around monitoring sites. In general, the intensity of noise is stronger than that of the MT signal, which can even be completely submerged by cultural noise, causing isolated points and distortions in apparent resistivity-phase curves, which are used to image geoelectrical structures. Therefore, it is necessary to improve the continuity and smoothness of the apparent resistivity-phase curves calculated based on impedance estimation. Two major types of methods are used in the impedance estimation of MT data according to the processing objects: frequency-domain methods and time-domain methods. One is an electric and magnetic field spectrum or power spectrum obtained by time–frequency transform, while the other is the time series of electric and magnetic fields and is directly observed.

Least squares estimation [4,5,6] and the robust method [7,8,9] are classic frequency-domain de-noising methods used in MT. There is a high correlation between horizontal MT signal components; thus, the least squares method, which uses a cross power spectrum, can eliminate the uncorrelated noise when there is no correlation between the signal and noise. However, the noise is always correlated with measured data; thus, the improvement of the traditional least squares method is biased. Based on the least squares method, the robust method was introduced for MT de-noising and can reduce the weight coefficients of noise through statistical analysis and coherency [7,10]. Smirnov (2003) proposed the criteria of the weight coefficient based on repeated median (RM), which is suitable for use in a strong-noise environment [11]. Chave (2014) introduced the maximum likelihood estimation method, which makes use of the distribution of noise data to calculate weights [12]. However, the robust method is insensitive to correlated noise and dependent on the quality of the data. In order to reduce noise weight coefficients, the amplitude of the signal must be far higher than that of the noise [13,14]. For correlated noise, Gamble et al. (1979) proposed the remote reference (RR) method. However, in industrialized areas such as ore mining areas, the quality of data is poor, meaning that any improvement made by the RR method is weak [15,16,17]. Oettinger et al. (2001) further introduced a signal-to-noise separation method with two remote reference sites to improve the de-noising ability of the RR method when used for a single site [18]. However, there has always been no solution to selecting the location of reference sites and the method is expensive.

In addition to the above-mentioned frequency-domain methods, time-domain methods can be used to improve the accuracy of impedance estimation. Empirical mode decomposition (EMD) [19,20] adaptively decomposes signals to intrinsic mode functions (IMFs) with different frequencies and a residual component. However, there the mode-mixing phenomenon occurs in EMD. The use of EEMD [21] and CEEMD [22] has been further proposed, which weakens the effect of mode-mixing, but there are also other problems, such as the existence of recovery error and the longer computing times needed. The mathematical morphological filtering method [23] and sparse decomposition algorithm [24] are effective only when the selection of parameters in the method is appropriate, which is affected by human factors that interfere greatly, and it is difficult to suppress multiple types of noise simultaneously. Morlet et al. (1982) introduced the method of continuous wavelet transform to signal analysis [25] and a series of mother wavelet transforms was further proposed [26,27,28,29]. Kumar and Foufoula-Georgiou (1997) proposed the application of discrete wavelet transform to geophysical signal processing [30]. Trad and Travassos (2000) introduced the wavelet threshold de-noising method based on the distribution of noise data [31]. Escalas et al. (2013) introduced polarization analysis for processing signals in wavelet domain in order to separate signals and noise [32]. However, the improvement of wavelet transform is also dependent on the selection of the parameters, such as the mother wavelet function and wavelet decomposition levels. Meanwhile, the most common way to select parameters for a method is through experience or repeated experiments, but this approach is influenced by subjective factors as well as being time-consuming.

Obviously, these methods are ineffective in improving the quality of data when the signal-to-noise ratio is low and the selection of the parameters used in methods is inappropriate. When the selected parameters are inappropriate, losses of signals can occur, or the methods may not even be capable of removing noise. Thus, we propose the use of a de-noising method based on discrete wavelet transform and iterative singular value decomposition (SVD). In this paper, the traditional SVD is improved by iterative loop that can extract noise more completely. Combined with discrete wavelet transform, the use of this method for the identification of various types of noise, including large-scale noise and impulse noise, is feasible. Meanwhile, multiscale dispersion entropy and phase space reconstruction are applied in this method to reduce the loss of signal and adaptively select parameters. To demonstrate the accuracy and stability of the proposed method, the proposed method is applied to process the MT field data of Linze area and Qilian area, China. Compared with traditional SVD decomposition, wavelet transform and the robust method, the proposed method can remove the various noise more thoroughly and obtain the useful MT response curves, which more truly reflect the subsurface electromagnetic structure.

2. Methods

2.1. Multiscale Dispersion Entropy

In general, compared with natural MT signals, cultural noise is expected to show a certain direction of polarization. Multiscale dispersion entropy (MDE) can identify the differences between the signal and cultural noise, and the MDE of cultural noise is smaller than that of signals. Compared with multiscale sample entropy (MSE), multiscale fuzzy entropy (MFE), and multiscale approximate entropy (MAE), the calculation of MDE is simpler and faster [33]. Multiscale dispersion entropy is defined as follows [34,35]:

(1): Time series u_j (j = 1, 2, …, N) are taken via a coarse graining preprocessing for non-overlapping segments x_b and are obtained as follows:

$x_{b} = \frac{1}{τ} \sum_{j = (b - 1) * τ + 1}^{b * τ} u_{j}, (1 \leq b \leq L)$

(1)

where τ is the scale factor, L = [N/τ] represents the length of the data segment, N is the length of time series, and [•] is the integral function.
(2): x_b is mapped to c classes from 1 to c. First, the normal cumulative distribution function (NCDF) is employed to map the segment to y_b from 0 to 1:

$y_{b} = \frac{1}{σ \sqrt{2 π}} \int_{- \infty}^{x_{b}} e^{\frac{- {(t - r m s)}^{2}}{2 σ^{2}}} d t$

(2)

where rms represents the root mean square and σ is the standard deviation (SD) of the segment, t is the length of observation time in each data segment. A linear algorithm z_b^c = round(c• y_b + 0.5) is used for y_b to z_b^c from 1 to c. Then, the embedding dimension m and time delay d are introduced to reconstruct z_i^m,c as follows:

$z_{i}^{m, c} = {z_{i}^{c}, z_{i + d}^{c} \dots z_{i + (m - 1) d}^{c}}$

(3)

where i = 1, 2…L $-$ (m $-$ 1)d. The number of possible dispersion patterns $π_{v_{0} \dots v}_{_{m - 1}}$ that can be assigned to each time-series z_i^m,c is equal to c^m, since the signal has m members and each member can be an integer from 1 to c.
Finally, for each c^m potential dispersion pattern, the relative frequency is defined as follows:

$p (π_{v_{0} \dots v}_{_{m - 1}}) = \frac{N u m b e r {i | i \leq L - (m - 1) d, \begin{matrix} z_{i}^{m, c} & h a s & t y p e & π_{v_{0} \dots v}_{_{m - 1}} \end{matrix}}}{L - (m - 1) d}$

(4)

where $p (π_{v_{0} \dots v}_{_{m - 1}})$ represents the number of dispersion patterns $π_{v_{0} \dots v}_{_{m - 1}}$ that are assigned to z_i^m,c divided by the total number of embedding signals with embedding dimension m. $v_{0} = z_{i}^{c}$ , $v_{1} = z_{i + d}^{c}$ …, $v_{m - 1} = z_{i + (m - 1) d}^{c}$ .
(3): The MDE is obtained as follows:

$M D E_{r m s} = E (x, τ, m, c, d) = [e_{1}, e_{2} \dots, e_{τ}]$

(5)

where $e (x, m, c, d) = - \sum_{π = 1}^{c^{m}} p (π_{v_{0} \dots v_{m - 1}}) \times l n (p (π_{v_{0} \dots v_{m - 1}}))$ represents the dispersion entropy based on Shannon’s entropy.

2.2. Phase Space Reconstruction

As a part of the chaos theory, phase space reconstruction is applied to indicate the number of real phase spaces in one-dimensional vectors, such as dynamical systems and time series, with differences seen between the correlation of data [36,37,38]. The algorithm is employed twice in the proposed method for (1) the selection of wavelet decomposition levels and (2) matrix construction before SVD decomposition.

Suppose U = {u₁, u₂ … u_N}, where N is the length of the signal, is reconstructed to different phase space vectors x_i = {u_i, u_i₊_τ … u_i _{+ (m−1)}_τ} (i = 1, 2, …, m) through the method of delay, where τ denotes time delay and m is the embedding dimension. Thus, the two parameters are necessary for the reconstruction. In general, various methods can be used to obtain τ, such as autocorrelation function [39], experience [40], and mutual information function [41]. However, compared with the above two methods, the third method is expected to process non-linear and non-stationary data. For the selection of m in this paper, false nearest neighbors (FNN) [42] are introduced to obtain m, where m must satisfy m > 2h + 1 (Takens theory) [43] and h represents the real dimension of attractors.

2.2.1. Mutual Information Function

Consider a general coupled system (S,Q), where [s, q] = [x(t), x(t + τ)] and the entropy of the system is H. When x is measured at time t, the average uncertainty of x at time t + τ is defined by averaging H(Q|s_i) over s_i as follows:

\begin{array}{l} H (Q | S) & = \sum_{i} P_{s} (s_{i}) H (Q | s_{i}) \\ = H (S, Q) - H (S) \end{array}

(6)

where H(Q|s) is the uncertainty of q given a measurement of s and H(S,Q) is the uncertainty of q in isolation. Thus, the measurement of s that reduces the uncertainty of q is:

\begin{array}{l} I (Q, S) & = H (Q) + H (S) - H (S, Q) \\ = I (S, Q) \end{array}

(7)

where I(S,Q) is the mutual information of coupled system (S,Q). Thus, mutual information can be generalized as follows:

I_{n} (X_{0}, X_{1}, \dots, X_{n}) = \sum_{j} (H (X_{j})) - H (X_{0}, X_{1}, \dots, X_{n})

(8)

when the one-dimensional vector is the time-delay reconstruction and the first minimum in I_n(τ) is employed as time delay τ.

2.2.2. False Nearest Neighbors

The nearest neighbor of x_i = {u_i, u_i₊_τ … u_i_+(m−1)_τ} by x_i⁰ = {u_i⁰, u_i₊_τ⁰ … u_i_+(m−1)_τ⁰} within a certain distance in m dimensions is denoted. Then, the square Euclidian distance is as follows:

R_{i} (m) = \sqrt{{\sum_{j = 0}^{m - 1} (u_{i + j} - u^{0}_{i + j})}^{2}}

(9)

When dimension m increases to dimension m + 1, the distance is expressed as:

R_{i} (m + 1) = \sqrt{\sum_{j = 0}^{m - 1} {(u_{i + j} - u^{0}_{i + j})}^{2} + {(u_{i + m} - u^{0}_{i + m})}^{2}}

(10)

Then, a₁(i, m) is given as follows:

a_{1} (i, m) = \frac{R_{i} (m + 1)}{R_{i} (m)}

(11)

For a₁(i, m) ≥ 10, the neighbor is identified as false; in contrast, m is the embedding dimension.

2.3. Discrete Wavelet Transform

Wavelet transform is defined as continuous wavelet transform (CWT) [44] and discrete wavelet transform (DWT) [30]. DWT has been introduced in the proposed method, in order to inversely transform the signal in the wavelet domain into a time series, which is available for further processing by SVD decomposition. Mother wavelet function and wavelet decomposition levels are the parameters used in wavelet transform. In MT signals, coifN is most suitable as the mother wavelet function, where N = 5 obtains a better result [45,46]. Meanwhile, the wavelet decomposition level in the proposed method is calculated using phase space reconstruction.

The definition of the DWT of u(t) is:

W_{f} (m, n) = λ_{0}^{- m / 2} \int u (t) ψ (λ_{0}^{- m} t - n t_{0}) d t

(12)

Ψ_{m, n} (t) = \frac{1}{\sqrt{λ_{0}^{m}}} Ψ (\frac{t - n t_{0} λ_{0}^{m}}{λ_{0}^{m}}) = λ_{0}^{- m / 2} Ψ (λ_{0}^{- m} t - n t)

(13)

Equation (13) represents the mother wavelet function. m and n are the integers. λ₀ is the scale. t is the temporal shift and t₀ is initial condition. W_f (m, n) denotes the wave coefficients, which satisfy the following expression:

A {‖ f ‖}^{2} \leq {\sum_{m} \sum_{n} | W_{f} (m, n) |}^{2} \leq B {‖ f ‖}^{2}

(14)

where A > 0, B < ∞ are both characteristic constants of the wavelet and the choices of λ₀ and t₀ are correlated with them [30]. Thus, the inverse wavelet transform is obtained as follows:

f (t) = \frac{2}{A + B} \sum_{m} \sum_{n} W_{f} (m, n) Ψ_{m, n} + γ

(15)

2.4. SVD Decomposition

Suppose that X = [x₁^T, x₂^T, … x_n^T]^T, where X∈R^m ^{× n}. The following formulation of the SVD method [47] is:

X = U S V^{T}

(16)

where

S = {\begin{matrix} {(d i a g (σ_{1}, σ_{2}, \dots σ_{P}), O)}^{T}, m > n \\ {(d i a g (σ_{1}, σ_{2}, \dots σ_{P}), O)}^{}, m < n \end{matrix}

, S∈R^m ^{× n} represents a diagonal matrix, O is the zero matrix, σ₁ ≥ σ₂ ≥ …≥ σ_p ≥ 0, and p = min(m, n), which are the singular values of X in descending order. U is an m × m matrix, of which the columns are orthonormal. V is an n × n matrix, of which the rows are orthonormal.

In general, there are two methods for transforming a one-dimensional vector to a matrix [48]: one involves non-overlapping the signal to N segments, while the other involves transforming the vector to a HERMIT matrix, where τ is 1. However, the methods above are special cases in the construction of the matrix and will not eliminate cultural noise effectively. Thus, the method of delay is applied in this paper; this selects the parameters τ and m with the consideration of data in different actual situations. Compared with artificially selected parameters, this method improves the precision and time-consuming characteristics. Meanwhile, we improve SVD decomposition with the iteration until the MDE of all components is more than the threshold, which is the termination criterion for the iteration, to achieve the complete separation of signal and noise.

σ^{i} = σ_{s}^{i} + σ_{n}^{i} (i = 1, 2, \dots, k)

(17)

Y = X_{s} + X_{n} = [\begin{matrix} U_{s} & U_{n} \end{matrix}] [\begin{matrix} S_{s} & 0 \\ 0 & S_{n} \end{matrix}] [\begin{matrix} V_{s}^{T} \\ V_{n}^{T} \end{matrix}]

(18)

where σ_sⁱ and X_s represent the effective signal section and σ_nⁱ and X_n are the noise section. The two-norm

f = {({‖ \sum \bar{X_{s}} ‖}_{2})}_{m i n}

is employed to estimate the signal and the reconstruction matrix of the signal is expressed as follows:

\bar{X_{s}} = U_{s} S_{s} V_{s}^{T}

(19)

where

\bar{X_{s}}

is the average of all the reconstruction matrices. Equations (15)–(18) are repeated until

M D E (\bar{X_{S}}) > ε

, where ε is the threshold, for which the value is the same as the MDE in Section 2.1, which is employed to filter noisy data.

The steps adopted in this paper are as follows (in Figure 1):

Overlap data to segments to ensure the continuity of the data;
Calculate the MDE for each data segment, for which the value below the threshold is the noisy data segment;
Apply phase space reconstruction to calculate the number of wavelet decomposition level in noisy data segments;
Perform discrete wavelet transform and discrete wavelet inverse transform for multiple components;
Decompose the components using iterative SVD to obtain the de-noised section;
Reconstruct the MT de-noising signal with the useful data segments given in Step 1 and the de-noising data segments given in Step 5, where the average value of the two segments is adopted for the overlapping data.

3. Synthetic Cases

To test the effectiveness of the proposed method, various simulated noise was randomly added into natural MT signals with a high signal-to-noise ratio, such as square wave noise, charge–discharge triangular wave noise, and impulse noise. In order to compare the recovery results of different methods, we introduced signal-to-noise ratio (SNR), mean square error (MSE), normalized cross coherence (NCC), and data error (E) to evaluate the de-noising signals [49,50]:

SNR = 20 l g \frac{{‖ y (n) ‖}_{2}}{{‖ y (n) - x (n) ‖}_{2}}

(20)

MSE = \sqrt{\frac{\sum_{n - 1}^{N} {(y (n) - x (n))}^{2}}{N}}

(21)

NCC = \frac{\sum_{n - 1}^{N} y (n) \cdot x (n)}{\sqrt{(\sum_{n - 1}^{N} y^{2} (n)) \cdot (\sum_{n - 1}^{N} x^{2} (n))}}

(22)

E = \frac{{‖ y (n) - x (n) ‖}_{2}}{{‖ y (n) ‖}_{2}}

(23)

S / N = 20 l g \frac{{‖ y (n) ‖}_{2}}{{‖ x (n) ‖}_{2}}

(24)

where y(n) represents the de-noised signal, x(n) denotes the original signal, ||•||₂ is a 2-norm function, and S/N denotes the level of added noise. The value of 0 dB indicates that the intensity of the added noise is much stronger than the useful signal, which means that the noise level is high. On the contrary, 40 dB represents that the intensity of added noise is low. When S/N is 40 dB, the intensity of the noise is only one percent of the useful signal, meaning that it has little influence on the useful signal.

As shown in Figure 2, various types of noise (S/N = 0) are added into the synthetic noisy signal. In Figure 2a, the natural MT signal is measured in Linze Province with a sampling rate of 15 Hz. Figure 2b shows the addition of square wave noise to the time series, Figure 2c denotes the addition of charge–discharge triangular wave noise, Figure 2d shows the addition of impulse noise, and Figure 2e shows the addition of the above various types of noise. All the noise is added in random positions. Obviously, the useful signal is completely submerged and its amplitude is far less than that of the noise. Figure 2f–j show the frequency spectrum of Figure 2a–e. The amplitude shown in Figure 2f is less than the amplitude shown in Figure 2g,h,j.

3.1. Entropy Analysis

To reduce the loss of signal, multiscale dispersion entropy (MDE) was used to preprocess the overall data, with noisy data being selected. The MDE was calculated for the useful signal and noisy signal shown in Table 1 in order to test whether the difference between different signals could be found. This indicated that there was a sharp gap between the MDE of noise-free data and that of noisy data, meaning that MDE is capable of distinguishing different signals. Furthermore, the reason why MDE is more suitable as an evaluation criterion of noisy data is also shown in Table 1. The multiscale approximate entropy (MAE) and multiscale fuzzy entropy (MFE) of different signals are approximate (the difference is less than 0.05); thus, it was difficult to differentiate the noise from the useful signal. However, there was an obvious difference between the useful signal and noisy signal in MDE and multiscale sample entropy (MSE). Figure 3 shows their stability in different conditions. The green dashed line shows the baseline of the MDE, the blue dashed line shows the baseline of the MSE, and both of the values are shown in Table 1. The sample points below the baseline are considered to be noisy data. The red solid line shows the MDE, while the black solid line shows the MSE. The MSE in the useful signal fluctuates sharply (in Figure 3a). Compared with the data points filtered by the MDE, there are fewer data points filtered by MSE (In Figure 3b–e). Therefore, the MDE is more suitable to be used as an evaluation method for distinguishing between signals and noise. Meanwhile, various types of noise have a different performance of MDE. The signal with square wave noise has the smallest MDE value, meaning that it is the easiest to be selected using MDE; the triangle wave is second and the MDE of the signal with the addition of impulse noise is the closest to that of the useful signal.

3.2. Parameter Calculation

After selecting noisy data, the parameters of the method should be calculated. In the proposed method, the wavelet decomposition level needs to be determined. The whole data set is divided into different data segments in order to calculate the MDE value for the evaluation of whether the data segment contains noise; thus, the number of noisy data segments is different for signals with the addition of different noise. There are six noisy data segments in the signal with square wave noise. There are 10 noisy data segments in other signals. The most common method used to determine the wavelet decomposition level is repeated experiments. In this method, the operator sets the value of the parameter one by one. To choose the best value, the four signal properties are calculated for each value. The result of the repeated tests is shown in Figure 4. When the number is 4 or 5, the improvement is obvious. Thus, the value is set as 4 or 5 by operators through repeated tests. To verify the reliability of the parameter calculated in the proposed method, the results determined by the proposed method in different data segments are shown in Table 2. These values are mainly around 4 or 5, and the same as those calculated by the common method. In other words, the values of the parameters calculated by the proposed method are accurate and this method is more efficient.

3.3. Performance Evaluation

In the section, we analyze the de-noising performance of the proposed method. Figure 5 shows that the results gained by extracting different types of noise using the proposed method. In Figure 5a–d, the signal is contaminated with different types of noise and the noisy data segments are completely selected by MDE with little loss of useful data around the noise. However, there is an inevitable loss of useful signal, which is overlapped in the noise, meaning that the noisy signal must be further processed by the proposed de-noising method. As shown in Figure 5e–h, the profiles of noise are extracted and the fluctuation of the profile of the noise is decreased after processing by the proposed method. Compared with processing all the data points, the proposed method with MDE can reduce the loss of useful signals. In Figure 5i–l, it is indicated that the frequency amplitude of the noise extracted by the proposed method is same as the amplitude of noisy data (see Figure 2g–j). Therefore, on the basis of reducing the useful signal, the proposed method is capable of removing both large-scale noise and impulse noise. However, compared with the triangular wave noise and impulse noise, the improvement of square wave noise achieved by the proposed method is better, with less loss of signals.

Figure 6a–j show the time series after eliminating various types of noise using the proposed method; the corresponding frequency spectrum is shown in Figure 6k–o. Due to the lack of noisy data, the result of Figure 6f is the same as that of the original signal (see Figure 6a). The time series after eliminating the square wave noise is shown in Figure 6g, where it can be seen that the profile of the useful signal is similar to that shown in Figure 6a. The results after suppressing the triangular wave noise and various other types of noise are shown in Figure 6h,j. It is evident that there are losses around 500 data point and 2900 data point (see red boxes and blue boxes), with fewer losses seen in Figure 6j. In other words, when there are various different types of noise in the MT data simultaneously, the improvement gained by removing triangular waves through the proposed method is better. Figure 6i shows the time series after eliminating the impulse noise, which is roughly the same as the noise-free data.

The improvement results after de-noising through singular value decomposition (SVD), wavelet transform (WT), and the proposed method (S-W) are shown in Figure 7. Compared with the results gained using the above two methods for different types of noise, it can be seen that the proposed method removes the profiles of noise more completely. For large-scale noise, SVD decomposition and wavelet transform find it difficult to suppress the interference, with residual noise being left (see blue lines and red lines in Figure 7a–c,j–l). For impulse noise, the proposed method also achieved better results than the other methods (In Figure 7d,e). However, compared with the time series seen after eliminating various types of noise in Figure 7j–l, it can be seen that only subtle differences are achieved from eliminating one type of noise through different methods. Thus, we introduced the use of a frequency spectrum to show the differences between the methods (see Figure 8). Figure 8d–f show the frequency spectrum after impulse noise was removed. The amplitude is similar because impulse noise has little influence on the frequency domain. For large-scale noise, obvious differences can be seen in the 0–1 Hz frequency band (see green lines in Figure 8a–c,g–i,j–l). The amplitude of the spectrum after SVD decomposition and wavelet transform were performed is greater than that after the proposed method was performed. In other words, it is difficult for the two above methods to completely suppress various types of noise.

As shown in Figure 9, four signal properties were calculated to evaluate the de-noising performance. The green line represents the proposed method, the red line shows the SVD decomposition, and the black line indicates the wavelet transform. With the increase in S/N (namely, the level of noise decreased), both the SNR and NCC increased correspondingly and the E and MSE gradually decreased. The values of the signal properties in the proposed method were always better than those of wavelet transform and SVD decomposition. Therefore, the proposed method has more advantages for improving the signal-to-noise ratio.

4. Implementation for MT Field Data

In this section, the MT filed data are applied to verify the practicability of the proposed method to improve the MT response, which is measured in Qilian County, Qinghai Province, with a sampling rate of 15 Hz and a distance between the two sites of 10 km. The survey area is quite far away from the urban area; thus, the data quality of most of the sites is good. However, site 380 is close to the industrialized area and the data are contaminated by various noises, resulting in a low signal-to-noise ratio.

We randomly selected from the sites with high-quality data, and then took site 300 and site 380 as examples. As shown in Figure 10, the results obtained after removing different noise levels are compared with the raw time series. In Figure 10a, it can be seen that the field data of site 380 are contaminated by square wave noise and triangular wave noise. Figure 10b shows the profile extracted by the proposed method, which indicates that the noise is completely abstracted and that the trend of the useful signal is recovered after processing using the proposed method (see Figure 10a,c). In Figure 10d, the signal-to-noise ratio of site 300 is high; the data are only contaminated by impulse noise, which can also be eliminated using the proposed method. Thus, the improvement obtained by the proposed method is effective, whether the data are measured with a high noise level or a low noise level.

Figure 11 shows the improvement of site 380 through SVD decomposition, wavelet transform, and the proposed method. However, when large-scale noise is subtracted in the previous two methods, there is still impulse noise remaining (see Figure 11b–d). Meanwhile, in the 0–2 Hz (Figure 10f–h) region, it can be seen that the amplitude of the SVD decomposition and wavelet transform is more than that of the proposed method, showing that the signal and noise separation in the proposed method is better and more effective.

The apparent resistivity phase curves of site 380, site 360, site 340, site 320, site 300, and site 280 processed using the robust method and the proposed methods are shown in Figure 12. There is are high-quality data for site 360, site 320, and site 300, which is less contaminated by noise. For these sites, the results of the robust method and the proposed method are almost the same, which indicates that the improvements made by the proposed method are reliable and reasonable. However, when the measured data are contaminated with a higher noise level (which is the case for site 340, site 280, and site 380), the robust method can only improve the fluctuation of some jump points, indicating that the improvement is weaker than that gained in the proposed method. For poor-quality data, the proposed method is able to recover the smoothness of the apparent resistivity phase curves. Although the improvement of the proposed method is more obvious for ρ_yx in site 340 and ρ_xy in site 280, the de-noising curves of two methods for ρ_xy in site 340 and ρ_yx in site 280 are overlapping, which further verifies the stability and effectiveness of the proposed method.

5. Conclusions

To solve the MT response distortions caused by high levels of noise, we propose the use of a de-noising method based on discrete wavelet transform and singular value decomposition (SVD) that consists of three sections. Firstly, multiscale dispersion entropy (MDE) selects noisy data to reduce the loss of useful signal. Then, phase space reconstruction calculates the values of parameters in the method. Finally, discrete wavelet transform and iterative SVD decomposition realize the signal–noise separation in MT data.

Compared with the traditional SVD decomposition and wavelet transform achieved in synthetic tests, the proposed method is more able to eliminate various types of noise without noise remaining. Meanwhile, the data measured in the Qilian area were used to test the reasonability of the proposed method. Our method improves the quality of the data measured in the Qilian area. Furthermore, based on the more continuous and smoother response curves gained using the proposed method than the robust estimation, we can deduce that the proposed method is more effective and more suitable for improving the quality of data with a low signal-to-noise ratio.

The proposed method has many advantages for processing poor-quality data; thus, it can be applied for the de-noising of MT data in environments with strong interference (such as industrialized areas). In comparison with conventional methods, the parameters of the proposed method are selected adaptively according to the characteristics of the data with no need for manual intervention, which reduces the amount of time needed.

Author Contributions

Conceptualization, R.Z. and J.H.; methodology, R.Z., J.H. and Z.G.; software, R.Z.; formal analysis, R.Z., J.H. and Z.G.; investigation, R.Z.; writing—original draft preparation, R.Z.; writing—review and editing, R.Z., J.H., T.L. and Z.G.; visualization, R.Z.; funding acquisition, J.H. All authors have read and agreed to the published version of the manuscript.

Funding

This paper is financially supported by the National Key Research and Development project(2017YFC0601305), the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (Grant No. 2019QZKK0701), the National Natural Science Foundation of China (41504076) and the Jilin Province Geological Exploration Fund Project (2018-19).

Conflicts of Interest

The authors declare no conflict of interest.

References

Oliver-Ocaño, F.M.; Gallardo, L.A.; Romo-Jones, J.M.; Perez-Flores, M.A. Structure of the Cerro Prieto Pull-apart basin from joint inversion of gravity, magnetic and magne-totelluric data. J. Appl. Geophys. 2019, 170, 103835. [Google Scholar] [CrossRef]
Xin, Z.; Han, J.; Liu, C.; Liu, L.; Tian, Y.; Yang, B.; Li, Z.; Mu, Q. Lithospheric structure and transformation of the Songnen and Jiamusi blocks in the eastern segment of the Central Asian Orogenic Belt revealed by 3D magnetotelluric inversion. Tectonophysics 2021, 816, 229015. [Google Scholar] [CrossRef]
Wu, Y.; Han, J.; Liu, Y.; Ma, G.; Han, F.; Yang, Y.; Liu, L.; Guo, L.; Guan, Y.; Zhang, Y.; et al. Metallogenic model of the Shuangjianzishan Ag-Pb-Zn district, Northeast China: Revealed from integrated geophysical investigation. Geosci. Front. 2021, 10132, (prepubilsh). [Google Scholar] [CrossRef]
Sims, W.E.; Bostick, F.X., Jr.; Smith, H.W. The estimation of magnetotelluric impedance tensor elements from measured data. Geophysics 1971, 36, 938–942. [Google Scholar] [CrossRef] [Green Version]
Gamble, T.D.; Goubau, W.M.; Clarke, J. Error analysis for remote reference magnetotellurics. Geophysics 1979, 44, 959–968. [Google Scholar] [CrossRef]
Jones, A.; Jödicke, H. Magnetotelluric transfer function estimation improvement by a coherence-based rejection technique. In SEG Technical Program Expanded Abstracts 1984; Society of Exploration Geophysicists: Houston, TX, USA, 1984; Volume 50, pp. 51–55. [Google Scholar]
Egbert, G.D.; Booker, J.R. The robust estimation of geomagnetic transfer functions. Geophys. J. Int. 1986, 87, 173–194. [Google Scholar] [CrossRef] [Green Version]
Chave, A.D.; Thomson, D.J. Bounded influence magnetotelluric response function estimation. Geophys. J. Int. 2004, 157, 988–1006. [Google Scholar] [CrossRef]
Sutarno, D.; Vozoff, K. Robust M-estimation of magnetotelluric impedance tensors. Explor. Geophys. 1989, 20, 383–398. [Google Scholar] [CrossRef]
Streich, R.; Becken, M.; Ritter, O. Robust processing of noisy land-based controlled-source electromagnetic data. Geophysics 2013, 78, E237–E247. [Google Scholar] [CrossRef]
Smirnov, M.Y. Magnetotelluric data processing with a robust statistical procedure having a high breakdown point. Geophys. J. Int. 2003, 152, 1–7. [Google Scholar] [CrossRef] [Green Version]
Chave, A.D. Magnetotelluric data, stable distributions and impropriety: An existential combination. Geophys. Geophys. J. Int. 2014, 198, 622–636. [Google Scholar] [CrossRef] [Green Version]
Weckmann, U.; Magunia, A.; Ritter, O. Effective noise separation for magnetotelluric single site data processing using a frequency domain selection scheme. Geophys. J. Int. 2005, 161, 635–652. [Google Scholar] [CrossRef] [Green Version]
Huang, N.E.; Shen, Z.; Long, S.R.; Wu, M.C.; Shih, H.H.; Zheng, Q.; Yen, N.C.; Tung, C.C.; Liu, H.H. The empirical mode decomposition method and the Hilbert spectrum for non-stationary time series analysis. Proc. Math. Phys. Eng. Sci. 1998, 454, 903–995. [Google Scholar] [CrossRef]
Clarke, J.; Gamble, T.D.; Goubau, W.M.; Koch, R.H.; Miracky, R.F. Remote-reference magnetotellurics: Equipment and procedures. Geophys. Prospect. 1983, 31, 149–170. [Google Scholar] [CrossRef] [Green Version]
Egbert, G.D. The robust multiple station magnetotelluric data processing. Geophys. J. Int. 1997, 130, 475–496. [Google Scholar] [CrossRef] [Green Version]
Ritter, O.; Junge, A.; Dawes, G.J. New equipment and processing for magnetotelluric remote reference observations. Geophys. J. Int. 1998, 132, 535–548. [Google Scholar] [CrossRef] [Green Version]
Oettinger, G.; Haak, V.; Larsen, J.C. Noise reduction in magnetotelluric time-series with a new signal-noise separation method and its application to a field experiment in the Saxonian Granulite Massif. Geophys. J. Int. 2001, 146, 659–669. [Google Scholar] [CrossRef] [Green Version]
Tang, J.T.; Hua, X.R.; Cao, Z.M.; Ren, Z.Y.; Duan, S.L. Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data. Chin. J. Geophys. 2008, 51, 603–610. (In Chinese) [Google Scholar]
Wang, Y.; Zhu, Y.; Wang, Q.; Yuan, S.; Tang, S.; Zheng, Z. Effective component extraction for hydraulic pump pressure signal based on fast empirical mode decomposition and relative entropy. AIP Adv. 2020, 10, 075103. [Google Scholar] [CrossRef]
Li, Y.; Gao, S.; Zhang, S.; He, H.; Xian, P.; Yuan, C. The baseline wander correction based on the improved ensemble empirical mode decomposition (EEMD) algorithm for grounded electrical source airborne transient electromagnetic signals. Geosci. Instrum. Methods Data Syst. 2020, 9, 443–450. [Google Scholar] [CrossRef]
Li, J.; Cai, J.; Tang, J.-T.; Li, G.; Zhang, X.; Xu, Z.-M. Magnetotelluric signal-noise separation method based on SVM–CEEMDWT. Appl. Geophys. 2019, 16, 160–170. [Google Scholar] [CrossRef]
Li, G.; Liu, X.; Tang, J.; Li, J.; Ren, Z.; Chen, C. De-noising low-frequency magnetotelluric data using mathematical morphology filtering and sparse representation. J. Appl. Geophys. 2020, 172, 103919. [Google Scholar] [CrossRef]
Li, J.; Liu, X.; Li, G.; Tang, J. Magnetotelluric noise suppression based on impulsive atoms and NPSO-OMP algorithm. Pure Appl. Geophys. 2020, 177, 5275–5297. [Google Scholar] [CrossRef]
Morlet, J.; Arens, G.; Fourgeau, E.; Giard, D. Wave propogation and sampling theory, 1, Complex signal and scattering in multi-layered media. Geophysics 1982, 47, 203–236. [Google Scholar] [CrossRef] [Green Version]
Daubechies, I. Orthonormal bases of compactly supported wavelets. Commun. Pure Appl. Math. 1988, 41, 901–996. [Google Scholar] [CrossRef] [Green Version]
Mallat, S.G. Multiresolution approximations and wavelet orthonormal bases of L2(R). Trans. Am. Math. Soc. 1989, 315, 69–87. [Google Scholar]
Meyer, Y. Wavelets and Operators; Cambridge University Press: Cambridge, UK, 1990. [Google Scholar]
Chui, C.K. An Introduction to Wavelets; Academic Press: Cambridge, MA, USA, 1995; Volume 2, pp. 50–61. [Google Scholar]
Kumar, P.; Foufoula-Georgiou, E. Wavelet analysis for geophysical applications. Rev. Geophys. 1997, 35, 385–412. [Google Scholar] [CrossRef] [Green Version]
Trad, D.O.; Travassos, M.J. Wavelet filtering of magnetotelluric data. Geophysics 2000, 65, 482–491. [Google Scholar] [CrossRef]
Escalas, M.; Queralt, P.; Ledo, J.; Marcuello, A. Polarisation analysis of magnetotelluric time series using a wavelet-based scheme: A method for detection and characterisation of cultural noise sources. Phys. Earth Planet. Inter. 2013, 218, 31–50. [Google Scholar] [CrossRef]
Liu, H.W. A Study on Feature Selection Algorithms Using Information Entropy. Ph.D. Thesis, Jilin University, Changchun, China, 2010. [Google Scholar]
Zhang, X.; Li, J.; Li, D.; Li, Y.; Liu, B.; Hu, Y. Separation of magnetotelluric signals based on refined composite multiscale dispersion entropy and orthogonal matching pursuit. Earth Planets Space 2021, 73, 1–18. [Google Scholar] [CrossRef]
Amezquita-Sanchez, J.P.; Mammone, N.; Morabito, F.C.; Adeli, H. A new dispersion entropy and fuzzy logic system methodology for automated classification of dementia stages using electroencephalograms. Clin. Neurol. Neurosurg. 2021, 201, 106446. [Google Scholar] [CrossRef]
Chen, K.; Han, B.T. A survey of state space reconstruction of chaotic time series analysis. Comput. Sci. 2005, 32, 67–70. (In Chinese) [Google Scholar]
Yoo, M.W.; Rontani, D.; Létang, J.; Petit-Watelot, S.; Devolder, T.; Sciamanna, M.; Bouzehouane, K.; Cros, V.; Kim, J.V. Pattern generation and symbolic dynamics in a nanocontact vortex oscillator. Nat. Commun. 2020, 11, 601. [Google Scholar] [CrossRef] [Green Version]
de Paula Viveiros, A.M. Non-orbital characterizations of strange attractors: Effective intervals and multifractality measures. Chaos Interdiscip. J. Nonlinear Sci. 2021, 31, 033139. [Google Scholar] [CrossRef] [PubMed]
Jiayu, L.; Yueke, W.; Zhiping, H.; Zhenkang, S. Selection of proper time-delay in phase space reconstruction of speech signals. Signal. Process. 1999, 15, 220–225. (In Chinese) [Google Scholar]
Broomhead, D.S.; King, G.P. Extracting qualitative dynamics from experimental data. Phys. Sect. D Nonlinear Phenom. 1986, 20, 217–236. [Google Scholar] [CrossRef]
Fraser, A.M.; Swinney, H.L. Independent coordinates for strange attractors from mutual information. Phys. Rev. A 1986, 33, 1134. [Google Scholar] [CrossRef] [PubMed]
Kennel, M.B.; Brown, R.; Abarbanel, H. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A 1992, 45, 3403–3411. [Google Scholar] [CrossRef] [Green Version]
Takens, F. Detecting strange attractors in turbulence, dynamical systems and turbulence. Lect. Notes Math. 1981, 898, 366–381. [Google Scholar]
Grossmann, A.; Morlet, J. Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM J. Math. Anal. 1984, 15, 723–736. [Google Scholar] [CrossRef]
Wei, D. Coiflet-Type Wavelets: Theory, Design and Applications. Ph.D. Thesis, University of Texas at Austin, Austin, TX, USA, 1998. [Google Scholar]
Carbonari, R.; D’Auria, L.; Di Maio, R.; Petrillo, Z. De-noising of magnetotelluric signals by polarization analysis in the discrete wavelet domain. Comput. Geosci. 2017, 100, 135–141. [Google Scholar] [CrossRef]
Kousika, N.; Premalatha, K. An improved privacy-preserving data mining technique using singular value decomposition with three-dimensional rotation data perturbation. J. Supercomput. 2021, 77, 10003–10011. [Google Scholar] [CrossRef]
Shao-Qing, Y.; Chuan-Ying, J. Two practical methods of phase space reconstruction. Acta Phys. Sin. 2002, 51, 2452–2459. (In Chinese) [Google Scholar]
Candès, E.J.; Wakin, M.B. An introduction to compressive sampling. IEEE Signal. Process. Mag. 2008, 25, 21–30. [Google Scholar] [CrossRef]
Li, G.; He, Z.; Tang, J.; Deng, J.; Liu, X.; Zhu, H. Dictionary learning and shift-invariant sparse coding de-noising for CSEM data combined with CEEMD. Geophysics 2021, 86, 1–52. [Google Scholar] [CrossRef]

Figure 1. Steps of the proposed processing scheme.

Figure 2. (a) Noise-free MT data, (b) data contaminated by square wave noise, (c) data contaminated by charge–discharge triangular wave noise, (d) data contaminated by impulse noise, (e) data contaminated by various noise. (f–j) is the frequency spectrum corresponding to (a–e). Noise-free MT data were collected in Linze Province with a sampling rate of 15 Hz. All the noise was randomly simulated.

Figure 3. Entropy value of noise-free data and noisy data at different data segments. (a) Noise-free MT data, (b) data contaminated by square noise, (c) data contaminated by charge–discharge triangular wave, (d) data contaminated by impulse noise, (e) data contaminated by various noise. The green dashed line is the baseline of MDE, and the blue dashed line is the baseline of MSE. The red solid line shows the MDE, and the black solid line shows the MSE.

Figure 4. The choice of the best wavelet decomposition level made through repeated tests. (a) shows SNR and MSE of the whole noisy signal at different wavelet decomposition levels, (b) shows NCC and E of the whole noisy signal at different wavelet decomposition levels. Red solid line indicates SNR, black solid line indicates MSE, red dashed line indicates NCC, and black dashed line indicates E.

Figure 5. Data detected by MDE and the noise time–frequency spectrum extracted by the proposed method. (a) shows the noisy data selected by MDE in the signal contaminated by square wave noise, (e) indicates the square wave noise extracted by the proposed method. (b,f) are the results of the signal contaminated by charge–discharge triangular wave. (c,g) denote the results of data contaminated by impulse noise. (d,h) represent the results of data contaminated by various noise types. (i–l) are the frequency spectrum corresponding to (e–h).

Figure 6. Time–frequency domain results after application of the proposed method. (a–j) shows the original and de-noised time series of noise-free MT data, data contaminated by square noise, data contaminated by charge–discharge triangular waves, data contaminated by impulse noise, data contaminated by various types of noise. (k–o) is the frequency spectrum corresponding to (f–j). Red boxes indicate the differences between the de-nosing signal and original signal at 500 data point. Blue boxes indicate the differences between the de-nosing signal and original signal at 2900 data point.

Figure 7. The de-noising results obtained after the use of different methods. Black lines represent the original signal as a reference, blue lines indicate the de-noising results of SVD, red lines denote the processing results of WT, purple lines show the results improved by the proposed method, completely covered by black lines. (a–c) show the time series of data contaminated by square noise after the use of SVD, WT, and the proposed method. (d–f) show the time series of data contaminated by impulse noise after the above different methods have been used. (g–i) show the time series of data contaminated by charge–discharge triangular waves after the above different methods have been used. (j–l) shows the time series of data contaminated by various types of noise after the above different methods have been used.

Figure 8. Frequency domain results after the use of the different methods. (a–c) show the spectrum of data contaminated by square noise after the use of SVD, WF, and the proposed method. (d–f) show the spectrum of data contaminated by impulse noise after the above different methods have been used. (g–i) show the spectrum of data contaminated by charge–discharge triangular waves after the above different methods have been used. (j–l) show the spectrum of data contaminated by various types of noise after the above different methods have been used.

Figure 9. (a) Recovery error (E), (b) NCC, (c) MSE, and (d) SNR of the signals recovered by different methods at different S/Ns. Black line indicates wavelet transform, red line denotes SVD decomposition, and green line shows the proposed method.

Figure 10. Time series segments of the real sites 380 and 360 in Qilian area with a sampling rate of 15 Hz. (a) Raw signal of Site 380, (b) noise extracted by the proposed method, and (c) signal de-noised by the proposed method. (d) Raw signal of site 360, (e) noise extracted by the proposed method, and (f) signal de-noised by the proposed method.

Figure 11. Time series segment of site 380 after the use of different methods. (a) Raw signal of site 380, (b) noise de-noised by SVD, (c) signal de-noised by WT, (d) noise de-noised by the proposed method, and (e–h) the frequency spectrum corresponding to (a–d).

Figure 12. Apparent resistivity phase curves of field sites. The red curves represent the results obtained using raw data. The green curves denote the results obtained using the data de-noised using the proposed method. The blue curves stand for the results obtained using the data de-noised by the robust method.

Table 1. Entropy values of noise-free data and noisy data.

Different Signals	MAE	MSE	MFE	MDE
Noise-free signal	0.2531	1.6493	0.1582	2.7568
With square wave	0.4519	0.9279	0.1369	1.1521
With triangle wave	0.3458	1.3142	0.1174	2.0091
With impulse noise	0.3379	1.5034	0.1041	2.3361
With various noise	0.4323	0.9013	0.1124	1.1055

Table 2. The wavelet decomposition level calculated by the proposed method for the different noisy data segments in the different signals (signals with different types of noise contain different numbers of noisy data segments).

Signals with Different Noise	Noisy Segment1	Noisy Segment2	Noisy Segment3	Noisy Segment4	Noisy Segment5	Noisy Segment6	Noisy Segment7	Noisy Segment8
Square wave	4	5	4	5	4	5	-	-
Triangle wave	5	5	4	5	4	4	4	4
Impulse noise	4	5	4	4	4	5	4	4
Various noise	4	5	4	4	4	4	4	4

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

© 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/).

Share and Cite

MDPI and ACS Style

Zhou, R.; Han, J.; Guo, Z.; Li, T. De-Noising of Magnetotelluric Signals by Discrete Wavelet Transform and SVD Decomposition. Remote Sens. 2021, 13, 4932. https://doi.org/10.3390/rs13234932

AMA Style

Zhou R, Han J, Guo Z, Li T. De-Noising of Magnetotelluric Signals by Discrete Wavelet Transform and SVD Decomposition. Remote Sensing. 2021; 13(23):4932. https://doi.org/10.3390/rs13234932

Chicago/Turabian Style

Zhou, Rui, Jiangtao Han, Zhenyu Guo, and Tonglin Li. 2021. "De-Noising of Magnetotelluric Signals by Discrete Wavelet Transform and SVD Decomposition" Remote Sensing 13, no. 23: 4932. https://doi.org/10.3390/rs13234932

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu