2.1 EMD method

The signal S(t) is decomposed into N-level IMF components and a residual component by the EMD method. The EMD involves the adaptive decomposition of signal S(t) by means of the sifting process. The term IMF is adopted because it represents the oscillation mode embedded in the data. The sifting process is defined by the following steps

1. Identify levels of decomposition N, and $r_{j-\mathrm{1}}\left(t\right)=S\left(t\right)$ as residual parameter.

2. Extract IMF_j:

   • a.

     all extrema of r_j−1(t);

   • b.

     interpolate local maxima and minima as the upper and lower envelopes separately by a cubic spline line, and compute “envelope” E_min(t) and E_max(t);

   • c.

     compute the average component $m\left(t\right)=(E_{\min }\left(t\right)+E_{\max }(t\left)\right)/\mathrm{2}$;

   • d.

     extract the detail component $D_i\left(t\right)=x\left(t\right)-m\left(t\right)$;

   • e.

     iterate steps (a) to step (d) on the detail component D(t) until the stopping criterion satisfies the threshold, dd<ε; once criterion is achieved, the D(t) is regarded as the effective IMF_j, which can also be regarded as the zero mean generally; calculate the stopping criterion:

     $$\begin{array}{}\text{(1)}& \mathrm{dd}=\sum {\displaystyle \frac{\mathrm{|}{D}_{i-\mathrm{1}}\left(t\right)-D_i\left(t\right){\mathrm{|}}^{\mathrm{2}}}{D_{i-\mathrm{1}}(t{)}^{\mathrm{2}}}}.\end{array}$$

3. Update the residual: $r_j\left(t\right)=r_{j-\mathrm{1}}\left(t\right)-{\mathrm{IMF}}_j\left(t\right)$; the residual is deemed as the input for a new round of iterations.

4. Repeat steps (2) and (3) until the value of j is equal to N; the stopping criterion threshold ε is set to between 0.2 and 0.3; the result of the sifting procedure is that S(t) will be decomposed into IMF_j(t), $j=\mathrm{1},\mathrm{\dots },N$ and residual r_N(t):

   $$\begin{array}{}\text{(2)}& S\left(t\right)=\sum _{j=\mathrm{1}}^N{\mathrm{IMF}}_j\left(t\right)+r_N\left(t\right).\end{array}$$

2.2 EEMD method

The EEMD method is an improved method based on the EMD method. Compared with the EMD method, the EEMD method resolves the mode mixing problem and achieves better performance by adding white noise to the original signal (Wu et al., 2009). For the EEMD method, the first step is to produce an ensemble of datasets by adding the finite amplitude σ of Gaussian distribution white noise to the original data. The σ stands for standard deviation of white noise. Then EMD method is applied to each realization of datasets to get IMF_i(t) with NE times repeatedly. The next step is that the expected ${\widehat{\mathrm{IMF}}}_j$ is obtained by averaging the respective components in each realization to compensate for the effect of the addition of Gaussian white noise.

where NE is the ensemble numbers. Finally, the result of the sifting procedure is that S(t) will be decomposed into ${\widehat{\mathrm{IMF}}}_j\left(t\right)$, $j=\mathrm{1},\mathrm{\dots },N$ and residual r_N(t).

where σ is set to between 0.05 and 0.2 and NE is set to between 100 to 400. In this paper, we set σ and NE to 0.1 and 200 respectively.

2.3 EEMD-AF method

The EEMD method is equivalent to a sifting filter which sifts out each IMF component from signal S(t) according to oscillations from fast to slow. The lower-index IMF component mainly contains fast oscillations; meanwhile the higher-index IMF component mainly contains slow oscillations. The baseline wander is expected to be captured by higher-index IMFs. The simple removal of several higher-index IMFs may introduce significant distortions in the reconstructed signal.

Thus, the baseline wander is distributed over the desired components in several higher-index IMFs. To suppress the baseline wander, this method introduces a group of adaptive low-pass filters to process the several higher-index IMFs successively. The sum of the output of these filters is regarded as the reconstructed baseline estimate. Finally, the de-noised signal can be obtained by subtracting an estimated baseline from the noisy signal.

First of all, we suppose the signal S(t) contained severe baseline wander. After processing by EEMD, S(t) will be decomposed into IMFs which be referred to as a_k(t).

where N is the number of IMFs. Then, it is important to find out how much IMFs contribute to the baseline wander. Denote this number value as M. The a_k(t) is processed from the higher to lower index by low-pass filter h_k(t). The output of filter is b_k(t).

where ∗ denotes the convolution. The h_k(t) is the Butterworth low-pass filter whose cut-off frequency is ω_k. As the IMF index decreases, fewer slow-oscillations component, but more signal components, are contained in each IMF. So, we design a group of adaptive low-pass filters whose cut-off frequency is decreased as the IMF index decreases. In other words, the first processing of filters is that the last IMF, a_N(t), convolves with the first low-pass filter whose cut-off frequency is ω_N(t). And the cut-off frequency decreases along with the k index filter decreases.

where α is set to between 0.1 and 0.99 and $k=N,\mathrm{\dots },\mathrm{2},\mathrm{1}$. By this means the filter output b_k(t) containing low-frequency components are extracted from each IMF. Because the algorithm has to be adaptive, the output can be used to determine the value of M as the condition of the reconstructed signal. According to the analysis of the procedure above, the amplitude of the baseline should gradually decrease as the IMF index decreases on filter output b_k(t). As a result, to determine the value of M, we consider using the evaluation coefficient function P_k as a stopping criterion where the SD(b_k) stands for standard deviation b_k.

where $k=\mathrm{1},\mathrm{2},\mathrm{\dots },N$. The operator flip is the flipped function that the data rearrange in the opposite direction. The evaluation coefficient threshold δ is set to between 0.05 and 0.1. If P_k<δ, the value $M=N+\mathrm{1}-k$. In this process, we set ω_N, α, and δ to 10, 0.9, and 0.01 respectively. The sum of the output of filtered IMFs whose index is from M+1 to N is regarded as the reconstructed baseline estimate.

Finally, to obtain a reconstructed de-noised signal, the baseline estimate is subtracted from the original signal.

3.1 Simulation data

In the GREATEM system, the transmitter injects a bipolar square wave current into the ground; meanwhile the receiver and front-mounted coil were installed on an aircraft in response to the vertical component of the induced electromotive force in a horizontal layered earth model (Nabighian et al., 1988). Response signals are related to the size and depth of the underground conductor, the length and current of transmitter, the equivalent area of receiving coil, the horizontal offset, the flight altitude, and so on. These parameters can be used to calculate the time domain response as a clean signal in the horizontal layer earth model for simulation. In Fig. 1, the model parameters as follows: the length of the transmitter line TX is 1000 m on the ground, the transmitter current is 10 A with 50 % duty cycles at 25 Hz, the equivalent area of receiving coil RX is 1000 m², the horizontal offset is 50 m, the flight altitude is 35 m, and the sample rate of the receiver is 32 kHz. In this paper, we consider a three-layer earth model where parameters are shown in Table 1. In the end, we calculated the corresponding time domain signal and the vertical response–decay curve on a three-layer earth model.

Figure 1GREATEM model based on a three-layer earth model. The TX is the transmitter line on the ground and the line length is 1000 m; the transmitter current is 10 A and the frequency is 25 Hz. The RX is the receiving coil and the equivalent area is 1000 m²; the offset is 50 m, the flight altitude is 35 m, and the sample rate of receiver is 32 kHz. The other model parameters are shown in Table 1.

Download

Table 1Parameters of the three-layer earth model.

Download Print Version | Download XLSX

Because of non-periodic and non-stationary characteristics of baseline wander, it is difficult to synthesize these noises from a simulation on the computer. The simulated signal is obtained by superimposing the clean signal on the field baseline wander measured by the inertial navigation system. Figure 2a is a simulated noisy signal which is obtained by adding baseline wander to a clean signal with the duration is 10 s. Figure 2b is the field baseline wander noise which is measured by the inertial navigation system with the duration is 10 s.

Figure 2The simulated noisy signal and baseline wander signal whose duration is 10 s. (a) The simulated noisy signals; (b) the field baseline wander measured by the inertial navigation system.

Download

3.2 Performance of the correction and analysis

In this paper, in order to quantitatively assess the de-noised quality between our method and other methods, we propose the signal-to-noise ratio (SNR) and mean-square error (MSE) to evaluate the correctional methods in Eqs. (11) and (12), in which S(n) is a synthetic clean signal, $\widehat{S\left(n\right)}$ is a de-noised signal, and L is the number of samples. The higher the SNR, the better the correctional effect; the lower the MSE, the better the fitting result.

There are comparisons of SNR and MSE obtained by three correctional methods for the noisy signal whose duration is 60 s. The correctional results are shown in Table 2. In the table, the term noisy signal means the simulated noisy signal before correction. The SNR value shows that the three methods have a remarkable improvement in signal quality. It is proved that EEMD-AF and the wavelet-based method had better correctional performance than EEMD. Quantitatively, the SNR of the EEMD-AF method is significantly close to the SNR achieved by the wavelet-based method. It is obvious that the EEMD-AF achieves a correction performance similar to the wavelet-based method.

Table 2Correctional result comparison of SNR of different methods.

Download Print Version | Download XLSX

For further analysis, the response–decay curve is related to the conductivity of underground geological bodies in the data processing of GREATEM. Besides the SNR comparison above, the original data are used for stacking, averaging, and extracting a secondary field to build a one-by-one test point along with the survey path. Then the number of time gates is 24 per test point, where the width of time gates increases approximately logarithmically.

To generate the anomaly curve profile image, we process the simulated noisy data and correctional data by stacking, averaging, extracting, and gating. The anomaly curve profile generated from the clean signal responses is shown in Fig. 3a, where the 24 parallel lines of time gates are represented along with the test point. And Fig. 3b is an anomaly curve profile generated from the noisy signal, where the fake anomaly from Gate 14 to 24 is identified clearly due to baseline wander affecting the horizontal layer model. From Gate 20 to 24, they are mixed with each other. After the processing using the wavelet-based method, the anomaly curve profile is shown in Fig. 3c, where the fake anomaly is not accurately represented and the curves are almost parallel to each other. Figure 3d is the anomaly curve profile using the EEMD-AF method; it is obvious that the paralleled curves between the gates are better than in the above method.

Figure 3Anomaly curve profile image generated from different processed datasets. The duration of raw data is 60 s and the stacking interval is 0.2 s; therefore the number of test points is 300. (a) The clean signal from the theoretical model; (b) the noisy signal containing baseline wander; (c) the correctional signal using the wavelet-based method; (d) the correctional signal using the EEMD-AF method. The label gate marked in each panel represents the number of time gates from 1 to 24. Every specific number of a time gate means a different time width, which increased logarithmically. The vertical coordinates are changed to logarithmic format in order to facilitate the display of detail of each gate in each panel.

Download

The comparison of SNR and MSE profiles produced by the datasets by different methods is illustrated in Fig. 4 along with the test point, where SNR and MSE of the noisy signal are marked as reference (black solid curve). In Fig. 4a, the black solid curve shows that the stacking and averaging may produce the improved SNR for a noisy signal. Quantitatively, the EEMD-AF and wavelet-based method yield an SNR which is significantly higher than the value achieved by the EEMD method. It is observed that there are fluctuations in the SNR using the wavelet-based method (red solid curve); meanwhile there are stabilities in the SNR using the EEMD-AF method (blue solid curve). And in Fig. 4b, the MSE curves indicate the same conclusion. Results from the comparison of the figures also show that the EEMD-AF method significantly outperforms the wavelet-based one for the suppression of non-stationary baseline wander.

Figure 4Comparison of SNR and MSE profiles produced by the datasets by different methods along with test point. (a) The contrast between SNR on different datasets; (b) the contrast between MSE on different datasets. The vertical coordinates are changed to logarithmic format in order to show a wide range of changes in MSE.

Download