# Time–frequency peak filtering for random noise attenuation of magnetic resonance sounding signal

Time–frequency peak filtering for random noise attenuation of magnetic resonance sounding signal Summary When measuring in a geomagnetic field, the method of magnetic resonance sounding (MRS) is often limited because of the notably low signal-to-noise ratio (SNR). Most current studies focus on discarding spiky noise and power-line harmonic noise cancellation. However, the effects of random noise should not be underestimated. The common method for random noise attenuation is stacking, but collecting multiple recordings merely to suppress random noise is time-consuming. Moreover, stacking is insufficient to suppress high-level random noise. Here, we propose the use of time–frequency peak filtering for random noise attenuation, which is performed after the traditional de-spiking and power-line harmonic removal method. By encoding the noisy signal with frequency modulation and estimating the instantaneous frequency using the peak of the time–frequency representation of the encoded signal, the desired MRS signal can be acquired from only one stack. The performance of the proposed method is tested on synthetic envelope signals and field data from different surveys. Good estimations of the signal parameters are obtained at different SNRs. Moreover, an attempt to use the proposed method to handle a single recording provides better results compared to 16 stacks. Our results suggest that the number of stacks can be appropriately reduced to shorten the measurement time and improve the measurement efficiency. Hydrogeophysics, Fourier analysis, Time-series analysis 1 INTRODUCTION Magnetic resonance sounding (MRS) is a geophysical method that can directly detect water because of its unique sensitivity to protons, which can be used to provide the depth profile of water content in the subsurface (e.g. Legchenko & Valla 2002; Hertrich 2008; Behroozmand et al.2015). MRS has been widely applied in hydrogeological investigations (Legchenko et al.2002). In recent years, this method has made rapid developments in logging (e.g. Walsh et al.2013), laboratory studies (e.g. Dlugosch et al.2016), and field measurements (e.g. Grunewald et al.2016). However, when detecting in the geomagnetic field (∼0.05 mT), the weak MRS signals are often corrupted by noise, which hinders the estimation of important parameters and restricts a more widespread use of the MRS method. The noise interferences that disturb the MRS signals mainly include spiky noise, power-line harmonic noise and random noise. The common de-noising workflow is as follows: (1) de-spiking, (2) harmonic noise cancellation and (3) random noise attenuation. To improve the signal-to-noise ratio (SNR) and better estimate the signal parameters, several effective noise-cancelling approaches are proposed, which are focused mainly on de-spiking and harmonic noise cancellation. Jiang et al. (2011) suggested a statistical stacking method to discard spikes based on the Romanovsky criterion. Dalgaard et al. (2012) implemented a spike detection algorithm based on a nonlinear energy operator. A spike from an electric fence was modelled in a study by Larsen (2016) as an impulsive excitation of a tuneable fourth-order bandpass filter and subtracted. Notch filtering (Legchenko & Valla 2003) was the most common and effective method for removing power-line harmonics in single-channel instruments. Then, based on the multichannel instrument presented by Walsh (2008), advanced techniques for noise cancellation were enabled and several methods were proposed (e.g. Walsh 2008; Dalgaard et al.2012; Müller-Petke & Costabel 2014). When there were multiple noise sources, a de-noising strategy that cancelled the power-line harmonics through a model-based approach combined with multichannel Wiener filtering was suggested by Larsen et al. (2014). It should be noted that both spiky noise from electric fences and power-line harmonic noise have specific mathematical expressions and can be methodically removed. Random noise commonly includes white noise, coloured noise, band-limited white noise, narrowband noise, etc. White noise is most common in electronic circuits with a constant power spectral density and homogeneous frequency components. In contrast, coloured noise with a wide band has a non-uniform power spectral density function (Gao 2011). Band-limited white noise and narrowband noise can be considered as white noise convolved with the impulse responses of an ideal low-pass filter and an ideal bandpass filter, respectively. The power spectral density of band-limited white noise is constant within a limited bandwidth of low frequency, and the power spectral density function of narrowband noise is limited to a narrow bandwidth. Random noise cannot be described by a specific mathematical expression, as it is difficult to remove. The methods involving random noise attenuation have rarely been reported in the current literature. At present, the most widely used method is stacking, that is, averages of multiple measurements (e.g. Legchenko & Valla 2002; Dalgaard et al.2012). However, the effect of stacking is limited in a noisy environment, and stacking is time-consuming. In the present study, we propose the use of time–frequency peak filtering (TFPF; Boashash 2003) for random noise attenuation in MRS. TFPF is a new method based on the theory of instantaneous frequency estimation and the property of noise accumulation. It was first developed by Arnold and Roessgen (Arnold & Boashash 1994; Arnold et al.1994; Roessgen & Boashash 1994) and was applied to newborn electroencephalograms (e.g. Boashash & Mesbah 2004) and seismic fields (e.g. Wu et al.2016). In the theoretical derivation part of this paper, the noisy signal is first transformed into the instantaneous frequency (IF) of an analytic signal, and the analysis of the analytic signal will then produce a time–frequency distribution and the underlying MRS signal is finally achieved by taking the peak of the time–frequency distribution. The theoretical considerations were validated using numerical simulations, and the signal parameters were well estimated at different SNRs. In addition, the method was tested on a synthetic MRS signal contaminated with real noise recordings, and the noise recordings were collected from two different field surveys. The test showed that the attempt to use the proposed method to handle a single recording provided a better result than stacking. Furthermore, the real field testing results were provided and discussed. We showed that random noise could be effectively suppressed, and an improved retrieval of the underlying MRS signal was acquired. 2 RANDOM NOISE IN MRS RECORDINGS The signals to be processed in this paper are the envelopes provided by MRS devices. The basic introduction to the envelope recordings is adequately described by Legchenko & Valla (2002). The original MRS signal is generated by the excited protons oscillating with the Larmor frequency, which is an exponentially decaying function of time as follows:   $$E(t,q) = {E_0}(q) \cdot {e^{{{ - t} / {T_2^*(q)}}}} \cdot \cos (2\pi {f_L}t + {\varphi _0}(q)),$$ (1)where t is the time sample, E0 is the initial amplitude of the MRS signal and is proportional to the water content, $$T_2^*$$ is the transverse relaxation time correlated to the pore geometry in an aquifer, fL is the Larmor frequency and may slightly vary from one place to another, φ0 is the phase shift between the recorded signal and the exciting current, q is the energizing pulse moment. The actual MRS signal is a multiexponential decay, but eq. (1) only explains a monoexponential, which can be considered an oversimplification. A high sampling frequency is required to measure the oscillating MRS signal, that is, dozens of kHz, so the amount of recorded data is large. To reduce the requirement for data acquisition and increase the measurement efficiency, a two-channel MRS signal detector followed by a low-pass filter is usually used to measure the envelope signal and the MRS devices, which record the envelope signal include JLMRS-I system (Jiang et al.2011), NumisPlus and NumisLite (Bernard 2006),   \begin{eqnarray} {E_{{\rm{Detection}}}}(t) &=& {E_0} \cdot {e^{{{ - t} / {T_2^*}}}} \cdot \cos (2\pi \Delta ft + {\varphi _0}) + {\varepsilon _X}(t)\nonumber\\ &&+\, i \cdot \left( {{E_0} \cdot {e^{{{ - t} / {T_2^*}}}} \cdot \sin (2\pi \Delta ft + {\varphi _0}) + {\varepsilon _Y}(t)} \right), \end{eqnarray} (2)where EDetection(t) denotes the envelope signal and Δf is the frequency offset of the Larmor frequency and the frequency of the excitation field. It is desirable that the frequency offset approach zero as much as possible. εX(t) and εY(t) are the additive measurement noises superimposed on the in phase and out phase of the MRS signal, which mainly consist of spiky noise, power-line harmonic noise and random noise. The random noise in the MRS originates from wide band background noise and receiver electronics noise. A field case of random noise attenuation by stacking is illustrated in Fig. 1. Figure 1. View largeDownload slide Field case: time-series and spectra of 1, 4, 8 and 16 stacks at one pulse moment, respectively. The left-hand column is the time-series. The blue lines display signals during the MRS experiments. The dashed black lines display noise recorded prior to the MRS signal measurement. The right-hand column is the spectra. The blue lines display the spectra of the signals, and the dashed black lines display the spectra of the noise that correspond to the left-hand column. Figure 1. View largeDownload slide Field case: time-series and spectra of 1, 4, 8 and 16 stacks at one pulse moment, respectively. The left-hand column is the time-series. The blue lines display signals during the MRS experiments. The dashed black lines display noise recorded prior to the MRS signal measurement. The right-hand column is the spectra. The blue lines display the spectra of the signals, and the dashed black lines display the spectra of the noise that correspond to the left-hand column. Fig. 1 shows the measurement results from 1, 4, 8, and 16 stacks at one pulse moment, which contained both signals and noise. The Larmor frequency was 2326 Hz. As widely used MRS working schemes, the field noise was directly recorded prior to the signal measurement. Here, spikes and power-line harmonics were removed prior to stacking. The time-series of the signal was corrupted by noise, and the frequency response of the signal was almost invisible from the spectral analysis for 1 stacked signal, as shown in Figs 1(a) and (b). Figs 1(c) and (d) illustrate that the noise was suppressed, and the signal was enhanced when the number of stacks was 4. When the number of stacks continued to increase, the noise level decreased. The comparisons between 8 stacks (row 3) and 16 stacks (row 4) showed that the SNR increased with the increasing stacks, but the random noise still remained. Assuming that si is the signal of the single recording and εi is the noise of the single recording, the SNR after N stacks is given as follows:   $${\rm {SNR}} = \frac{{\sum\limits_{i = 1}^{{N_s}} {{s_i}} }}{{\sqrt {\sum\limits_{i = 1}^{{N_s}} {{\varepsilon _i}^2} } }} = \frac{{{N_s} \cdot s}}{{\sqrt {{N_s}} \cdot \varepsilon }} = \sqrt {{N_s}} \cdot \frac{s}{\varepsilon },$$ (3)where $$\sum\limits_{i = 1}^{{N_s}} {{s_i}}$$ and $$\sqrt {\sum\limits_{i = 1}^{{N_s}} {{\varepsilon _i}^2} }$$ are the sum of the signals and noise intensity for Ns stacks, respectively; s and ε are the mean magnitude of the signals and noise for Ns stacks, respectively; $$\frac{s}{\varepsilon }$$ is the SNR before stacking, so that the SNR increases $$\sqrt {{N_s}}$$ times (where Ns is the number of multiple measurements). This relationship is obtained and shown in Fig. 2. Figure 2. View largeDownload slide Schematic of the signal-to-noise ratio with increasing number of stacks. Figure 2. View largeDownload slide Schematic of the signal-to-noise ratio with increasing number of stacks. Fig. 2 shows the schematic of the SNR with the increased number of stacks. This observation implies that the SNR cannot always be increased rapidly with the number of stacks and it has diminishing returns. From 1 stack to approximately 32 stacks, the SNR notably and rapidly increases; from 32 stacks to 100 stacks, it steadily increases. Then, the SNR increases very slowly after 100 stacks; that is, increasing the stack numbers is meaningless because the measurement takes a very long time. Here, we take 100 stacks as an example. If the multipulse Carr–Purcell–Meiboom–Gill (CPMG) sequence is used in the MRS measurement (Grunewald & Walsh 2013) and a full set of phase cycles is composed of four individual recordings, a complete CPMG data set takes 20 s. To obtain the accurate hydrologic information, 16 pulse moments need to be carried out. If each full cycle were repeated 100 times for noise reduction, the overall measurement time would be 12 or 13 h, including the time required to move the equipment and place the measurement coils in the field. Therefore, collecting multiple recordings to obtain the best SNR is time-consuming, and the method to average the time-series is not sufficient when the proton magnetic resonance response is covered by random noise. A more effective method is required. Thus, we developed the TFPF method, which handles a single recording to replace the stacking method for random noise attenuation. 3 METHODOLOGY 3.1 Time–frequency peak filtering algorithm Implementation of the TFPF algorithm is simple, whereby the noisy signal to be filtered is first encoded as the IF of a frequency-modulated (FM) analytic signal, and the analytic signal plays a key role in the TFPF algorithm. Analysing the analytic signal with the pseudo-Wigner-Ville distribution (PWVD) will produce a time–frequency distribution exhibiting significant energy concentration around the IF of the analytic signal. The IF is then estimated by taking the peak of the time–frequency distribution to recover the desired signal without losing useful information. When implementing the TFPF algorithm with PWVD, two special steps are introduced to improve the performance, including endpoint extension and scaling. The practical application of the TFPF algorithm is illustrated in Fig. 3, and a detailed description of the TFPF scheme is presented in the following steps. Figure 3. View largeDownload slide Flowchart of the TFPF algorithm. The TFPF algorithm consists of two stages. In the first stage, the noisy signal is encoded as the instantaneous frequency of the analytic signal after applying endpoint extension and scaling. In the second stage, the time–frequency distribution of the analytic signal is calculated using the PWVD, then take the peak of the PWVD of the analytic signal, followed by IF conversion, inverse scaling and endpoint rejecting, the filtered signal without random noise interference is finally obtained. Figure 3. View largeDownload slide Flowchart of the TFPF algorithm. The TFPF algorithm consists of two stages. In the first stage, the noisy signal is encoded as the instantaneous frequency of the analytic signal after applying endpoint extension and scaling. In the second stage, the time–frequency distribution of the analytic signal is calculated using the PWVD, then take the peak of the PWVD of the analytic signal, followed by IF conversion, inverse scaling and endpoint rejecting, the filtered signal without random noise interference is finally obtained. Let x(m) be the random-noise-corrupted MRS signal, which is modelled by the following equation:   $$x(m) = s(m) + \varepsilon (m),$$ (4)where m is the discrete time sample, m = 1, 2, …, M, M is the length of x(m), s(m) is the desired MRS signal and ε(m) is the random noise. 3.1.1 First stage: encoding First step: extend the endpoints of the sampled signal. To obtain the unbiased estimation of the desired signal, an appropriate window length is selected during PWVD. However, the available sampling points at the boundaries are less than the window length, and the energy would not concentrate at the boundaries of the analytic signal's IF on the time–frequency plane. As a result, the values of the estimated signal at both ends are deviated from the real values when the observed dataset is directly processed. To avoid endpoint distortion, endpoint extension should be performed, and the extended samples are required to be larger than half of the window length. Here, we chose a simple way to extend the observed dataset with a length of p samples at both ends, which is given as follows:   $$x'(n) = \left\{ \begin{array}{@{}l@{}} x(1),1 \le n \le p\\ x(n - p),p + 1 \le n \le p + M\\ x(M),p + M + 1 \le n \le 2p + M \end{array} \right.,$$ (5)where x'(n) is the signal after endpoint extension and 1 ≤ n ≤ 2p + M. Assuming that the window length W = 2L + 1, p ≥ L + 1 is required. Second step: scale the amplitude of the extended signal. Modulating the unscaled signal will lead to aliasing, which produces discontinuities in the estimated IF at the frequency boundaries of the time–frequency plane. To process all signals with different amplitudes, scaling should be performed before encoding to provide suitable frequency limits on the encoded signal. For easy understanding, the signal which is then encoded as the IF of the analytic signal is sampled at a normalized frequency of 1 Hz, the amplitude of the signal should range from half of the normalized frequency (0.5) to 0:   $${x_c}(n) = (a - b)\frac{{x'(n) - \min [x'(n)]}}{{\max [x'(n)] - \min [x'(n)]}} + b,$$ (6)where xc(n) is the scaled signal. After applying amplitude scaling, parameters a and b are the maximum and minimum of the discrete noisy signal, respectively, which satisfies the constraint 0 ≤ b < a ≤ 0.5. To preserve the signal details to make the amplitude of the scaled signal vary as large as possible, we normally let a = 0.5 and b = 0. Third step: encode the scaled signal as the IF of an analytic signal. According to the notion of IF (e.g. Goswami & Hoefel 2004), the IF of a given analytic signal z(t) = ej2πϕ(t) is as follows:   $${f_z}(t) = \frac{1}{{2\pi }}\frac{{{\rm{d}}\phi (t)}}{{{\rm{d}}t}},$$ (7)where ϕ(t) is the instantaneous phase and fz(t) is the IF of the analytic signal, and thus, the analytic signal can also be expressed as follows:   $$z(t) = {{\rm{e}}^{j2\pi \int_{{ - \infty }}^{t}{{{f_z}(\lambda ){\rm{d}}\lambda }}}}.$$ (8) Likewise, the discrete scaled signal xc(n) is encoded as the IF of the unit-amplitude FM analytic signal z(n), given as follows:   $$z(n) = {{\rm{e}}^{j2\pi \sum\limits_{\lambda = 0}^n {{x_c}(\lambda )} }}.$$ (9) 3.1.2 Second stage: IF estimation First step: complete the PWVD of the analytic signal. Calculate the time–frequency distribution of the analytic signal using PWVD, the discrete PWVD of the analytic signal z(n) can be expressed as follows:   $${\rm{PWV}}{{\rm{D}}_z}(n,k) = 2\sum\limits_{l = - L}^L {{\rm w}(l)z(n + l){z^*}(n - l){{\rm{e}}^{ - j4\pi kl}}} ,$$ (10)where w(l) is a window function, the symbol * denotes the complex conjugate, k is the discrete frequency sample, and 2L + 1 is the window length (W) considered in the PWVD discrete time implementation. It is important that the signal be linear to obtain an unbiased estimation (Boashash & Mesbah 2004). However, the MRS signal is a nonlinear function of time, so an appropriate window length should be selected to make the signal as close to linear as possible across the window length. The relationship related to the TFPF window length is given for the case of signal estimation, which is expressed as follows:   $$W \le \frac{{0.634{f_s}}}{{\pi f}}.$$ (11) Eq.(11) gives the maximum window length as a function of sampling rate (fs) and signal frequency (f). For the signals generated by protons, the frequency of the envelope signal would not exceed the range of 5 to −5 Hz (Strehl 2006). To reduce the bias, we select a smaller window length, and the signal frequency of 5 Hz is always chosen to calculate the window length. In addition, the sampling rate is normally four times the Larmor frequency. Without a loss of signal characteristics, resampling of the time-series is always needed, for example, the resampling rate of the envelope signal is a quarter of the Larmor frequency (Wang 2009). The random noise that corrupted the envelope signals can be attenuated effectively by the TFPF method, but it cannot process the oscillating signals well at present because the oscillating NMR signals recorded at dozens of kHz sampling rates have a high nonlinearity and the IF estimations are biased. Signal preservation and noise suppression are difficult to weigh. If the selected window length is large, the signal component will be lost. Conversely, the noise will not be better suppressed. This is the limitation of the present method. Second step: take the frequency samples of the IF estimation. Take the peak of PWVDz(n, k) to estimate the IF as follows:   $${\hat{k}_z}(n) = {\arg _k}\max [{\rm{PWV}}{{\rm{D}}_z}(n,k)],$$ (12)where $${\hat{k}_z}(n)$$ is the frequency sample of the IF estimation, which can be considered a recovery of the desired signal based on the standard time–frequency peak-detection algorithm (Boashash 1992). The subscript k for $${\arg _k}$$ means a discrete frequency sample; when maximizing the PWVD over the frequency, the peak of PWVD is found for each value of k. Third step: convert the frequency samples to the IF estimation. Because the values of $${\hat{k}_z}(n)$$ are composed of the frequency samples, we should convert the estimated frequency samples as follows:   $${\hat{s}_c}(n) = ({\hat{k}_z}(n) - 1)\frac{{{f_s}}}{N},$$ (13)where $${\hat{s}_c}(n)$$ is the estimated IF after conversion; $$\frac{{{f_s}}}{N}$$ is the frequency resolution; fs is the sampling rate and N is the number of frequency bins during the discrete PWVD. Fourth step: inverse scaling of the signal after IF conversion. The signal after IF conversation ranges from a to b, the order of the desired signal is recovered by an inverse scaling operation:   \begin{eqnarray} \hat{s}'(n) = \frac{{({{\hat{s}}_c}(n) - b)(\max [x'(n)] - \min [x'(n)])}}{{a - b}} + \min [x'(n)],\nonumber\\ \end{eqnarray} (14)where $$\hat{s}'(n)$$ is the desired signal with a real amplitude. Fifth step: reject the extended endpoints. The final filtered result $$\hat{s}(m)$$ shall be obtained by rejecting the extended endpoints given as follows:   $$\hat{s}(m) = \hat{s}'(m + p),1 \le m \le M.$$ (15) 3.2 A complete signal example In what follows, we show a complete signal example to explain each step of the TFPF presented in the previous section. A synthetic MRS envelope signal with the parameters E0 = 200 nV and T2* = 150 ms was generated and had been corrupted by random noise with a standard deviation of 50 nV (the mean value of the artificial random noise is 0 nV in all the following experiments), making the SNR 7.11 dB. Then, the TFPF method was performed to mitigate the random noise, and the basic procedure is illustrated in Fig. 4. Figure 4. View largeDownload slide Example of the basic procedure for the random noise attenuation. (a) The original noisy signal. (b) Endpoint extension. (c) Amplitude scaling. (d) Time–frequency implementation. (e) IF estimation. (f) IF conversion. (g) Inverse scaling. (h) Endpoints rejecting. Figure 4. View largeDownload slide Example of the basic procedure for the random noise attenuation. (a) The original noisy signal. (b) Endpoint extension. (c) Amplitude scaling. (d) Time–frequency implementation. (e) IF estimation. (f) IF conversion. (g) Inverse scaling. (h) Endpoints rejecting. The random-noise-corrupted MRS signal is shown in Fig. 4(a). The length of the observed signal was 145 (Larmor frequency fL = 2320 Hz, the resampling rate was set as a quarter of the Larmor frequency, and the recording time was 250 ms). In the first stage, we dealt with the noisy signal by extending the endpoints according to eq. (5), a length of eight samples were extended at both ends of the original noisy signal, which was shown in Fig. 4(b) (first step). Then, signal scaling was performed before encoding according to eq. (6). Assuming that the extended signal, which was regarded as the IF of the analytic signal was sampled at a normalized frequency of 1 Hz, Fig. 4(c) shows the amplitude of the scaled signal was between 0.5 and 0, where we selected the scaling parameters a = 0.5 and b = 0 (second step). Finally, the scaled signal was encoded as the IF of a unit-amplitude FM analytic signal according to eq. (9). Thus, the analytic signal was obtained, and the analytic signal was important in the TFPF scheme (third step). In the second stage, the PWVD was first performed on the analytic signal according to eq. (10). The window length (W) used in the discrete PWVD was 15. A time–frequency distribution exhibiting significant energy concentration around the IF of the analytic signal is produced, which resembled the desired MRS signal as shown in Fig. 4(d). However, the random noise was scattered over the time–frequency plane due to its stochastic properties, and this yielded an appropriate separation between signal and noise components, as well as a good estimation without the effect of random noise (first step). Subsequently, the frequency samples were obtained by maximizing the PWVD along the frequency axis according to eq. (12). The number of frequency bins (N) we selected during PWVD was 128, so the maximum value of y-axis was 128 as shown in Fig. 4(e) (second step). Then, the estimated frequency samples were multiplied by the frequency resolution (fs/N) according to eq. (13), the converted signal (Fig. 4f) was generated with the amplitude ranged from 0.5 to 0 (third step). Subsequently, inverse scaling was conducted according to eq. (14), the real amplitude of the signal was recovered (Fig. 4g, fourth step), followed by the endpoints being rejected according to eq. (15) (fifth step), the filtered signal was finally obtained with parameters E0 = 201.8 ± 2.1 nV, T2* = 149.97 ± 1.60 ms and the SNR increased to 28.02 dB, which is shown in Fig. 4(h). 4 NUMERICAL SIMULATIONS In this section, to validate the proposed method, numerical experiments were performed to model the monoexponential decaying signals in the artificial random noise and real noise recordings. To quantify the effects of the random noise attenuation, a fitting algorithm described in Legchenko & Valla (1998) was used to estimate the signal parameters. The performance of the proposed method was also evaluated based on the SNR and mean absolute percentage error (MAPE). Here, as a scale-dependent error metric (Hyndman & Koehler 2006), the MAPE was quoted to measure the quality of preserved information in the de-noised signal. The MAPE is defined as follows:   $$MAPE\,{\rm{ = }}\,\frac{1}{M}\sum\limits_{m = 1}^M {\left| {\frac{{\hat{s}(m) - s(m)}}{{s(m)}} \times 100} \right|} ,$$ (16)where M is the length of samples, s(m) is the ideal signal and $$\hat{s}(m)$$ is the de-noised signal. 4.1 Synthetic example for artificial noise A monoexponential decaying signal with the parameters E0 = 200 nV, T2* = 150 ms, Δf = 0 Hz and φ0 = 1.05 rad were generated. To make the synthetic data comparable with field conditions, the added noise was simulated as follows (Behroozmand et al.2013):   \begin{eqnarray} x(m) = s(m) + {\varepsilon _{G(0,1)}}(m) \cdot {\left[\varepsilon _{{\rm{uni}}}^2(m) + {\left(\frac{{{\varepsilon _{{\rm{bac}}}}(m)}}{{s(m)}}\right)^2} \right]^{\frac{1}{2}}} \cdot s(m),\nonumber\\ \end{eqnarray} (17)where m is the discrete time sample, x(m) and s(m) denote the noisy signal and ideal signal, respectively. εG(0, 1)(m) is the Gaussian noise with a mean value of 0 and a standard deviation of 1; εuni(m) is uniform noise, which can be considered non-specified noise contributions such as structural noise; and εbac(m) is the background noise. Fig. 5 shows the results of the TFPF method for random noise attenuation under different noise conditions. To simulate MRS signal 1, the monoexponential decaying signal was corrupted by a Gaussian noise distribution with a standard deviation of 30 nV and a uniform relative noise of 3 per cent of the data values, which were assigned to Vbac and STDuni in eq. (17). The SNR was 10.81 dB. Then, the TFPF method was employed to process the noisy signal. Panel (c) showed the PWVD of the analytic signal, the observation implied that a significant energy concentration was produced, and the signal could be estimated by maximizing the PWVD of the analytic signal along the frequency axis. The signal, regarded as the IF of the analytic signal, was sampled at a normalized frequency of 1 Hz, and the frequency axis of the time–frequency plane ranged from 1 to 0 Hz. Figs 5(a) and (b) illustrated the processing results in the time and frequency domains, respectively. The random noise component was eliminated, and the signal component was preserved. The estimated parameters, SNR and MAPE after TFPF, are reported in Table 1. We observed that with TFPF, the SNR increased to 28.32 dB, the estimated E0 was 199.3 ± 2.1 nV, T2* was 151.2 ± 1.5 ms, Δf was −0.01 ± 0.02 Hz, φ0 was 1.09 rad, and MAPE was 4.8 ± 1.4 per cent, respectively. To test the effectiveness of this method for a high noise level, the Gaussian noise was increased to a standard deviation of 150 nV, and the uniform relative noise was 5 per cent of the data values, which made the SNR of the signal 2 −3.16 dB. The middle row showed the outcome of the TFPF method, where the SNR increased to 26.27 dB. The estimated MRS signal parameters E0, T2*, Δf and φ0 were 201.2 ± 2.8 nV, 149.99 ± 2.05 ms, −0.03 ± 0.05 Hz and 1.11 rad, respectively, and the calculated MAPE was 5.7 ± 1.8 per cent. Furthermore, a serious noise disturbance was simulated as signal 3. The ideal signal was contaminated with Gaussian noise, the standard deviation was 300 nV with a uniform relative noise of 10 per cent of the data values, and the SNR was −10.66 dB. We observed that the signal was submerged in a high noise level from the time-series in Fig. 5(g). The frequency offset of 0 Hz was hardly observed and indistinctly emerged from the spectrum, as shown in Fig. 5(h). The fidelity of the MRS signal was severely degraded by random noise. However, the de-noised signal showed that the TFPF had an excellent performance because the high-level random noise did not produce an energy concentration in the time–frequency plane, as shown in Fig. 5(i). The high-level random noise was mitigated without losing the signal information, and the SNR increased to 21.73 dB; the estimated E0 was 198.9 ± 4.3 nV, T2* was 146.7 ± 5.1 ms, Δf was 0.04 ± 0.07 Hz, φ0 was 0.94 rad, and MAPE was 9.2 ± 3.3 per cent. In addition, we increased the noise level and tested this method. The results showed that the underlying signals could be clean recovered in the noise level down to an SNR of −17 dB. Figure 5. View largeDownload slide Numerical experiments of modelling the monoexponential synthetic signal added to three different artificial random noise recordings. The simulated signals with three different SNRs (signal 1 = 10.81 dB, signal 2 = −3.16 dB and signal 3 = −10.66 dB) are shown in the upper, middle and lower rows, respectively. The time-series and spectra of the envelope signals are displayed in the left-hand and middle columns, respectively. The right-hand column shows the PWVD of the FM analytic signal. Figure 5. View largeDownload slide Numerical experiments of modelling the monoexponential synthetic signal added to three different artificial random noise recordings. The simulated signals with three different SNRs (signal 1 = 10.81 dB, signal 2 = −3.16 dB and signal 3 = −10.66 dB) are shown in the upper, middle and lower rows, respectively. The time-series and spectra of the envelope signals are displayed in the left-hand and middle columns, respectively. The right-hand column shows the PWVD of the FM analytic signal. Table 1. Estimated parameters, SNR and MAPE performance of TFPF on synthetic MRS signals superimposed in different artificial random noise recordings. The synthetic MRS signal has E0 = 200 nV, T2* = 150 ms, Δf = 0 Hz and φ0 = 1.05 rad.   E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  Signal 1 (SNR = 10.81 dB)  199.3 ± 2.1  151.2 ± 1.5  −0.01 ± 0.02  1.09  28.32  4.8 ± 1.4  Signal 2 (SNR = -3.16 dB)  201.2 ± 2.8  149.99 ± 2.05  −0.03 ± 0.05  1.11  26.27  5.7 ± 1.8  Signal 3 (SNR = -10.66 dB)  198.9 ± 4.3  146.7 ± 5.1  0.04 ± 0.07  0.94  21.73  9.2 ± 3.3    E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  Signal 1 (SNR = 10.81 dB)  199.3 ± 2.1  151.2 ± 1.5  −0.01 ± 0.02  1.09  28.32  4.8 ± 1.4  Signal 2 (SNR = -3.16 dB)  201.2 ± 2.8  149.99 ± 2.05  −0.03 ± 0.05  1.11  26.27  5.7 ± 1.8  Signal 3 (SNR = -10.66 dB)  198.9 ± 4.3  146.7 ± 5.1  0.04 ± 0.07  0.94  21.73  9.2 ± 3.3  View Large 4.2 Synthetic example for real noise In this subsection, another synthetic example was shown. A monoexponential decaying signal with the parameters E0 = 75 nV, T2* = 100 ms, Δf = 0 Hz and φ0 = −0.53 rad, given by eq. (2), were added to noise-only recordings from two different sites. The noise levels at site 1 were smaller than those at site 2; 16 stacks at one pulse moment for site 1 and 32 stacks at one pulse moment for site 2 were performed. We conducted two MRS experiments at each site and collected the noise-only recordings. The actual noise interferences included spiky noise, power-line harmonic noise and random noise. The serious spiky noise was found in almost all recordings. Hence, a de-spiking algorithm and a harmonic noise cancellation algorithm (Jiang et al.2011) were first used in these original noisy recordings. After de-spiking and harmonic noise cancellation, the random noise remained. Then, we arbitrarily selected a single recording from the processed multiple recordings and applied the TFPF method to suppress the residual random noise. We compared this result with stacking. Fig. 6 shows a comparison of the traditional stacking method and TFPF method for random noise attenuation. The panels in the left-hand column illustrated that the random noise component still remained after stacking, whereas good results were obtained by the TFPF method, which handled the single recordings. The processed results in the middle column showed that the residual random noise was suppressed, and the decaying signal curve appeared close to the ideal signal. Compared with the traditional stacking results, the proposed method enabled us to recover the desired MRS signal even at a low SNR. Figure 6. View largeDownload slide Comparison of the stacking and the TFPF method for the residual random noise attenuation. The panels in the left-hand column display the de-noising results by stacking. The middle column shows the results before and after applying the TFPF method for an arbitrarily selected single recording. The right-hand column shows the spectra of the filtered signals. Figure 6. View largeDownload slide Comparison of the stacking and the TFPF method for the residual random noise attenuation. The panels in the left-hand column display the de-noising results by stacking. The middle column shows the results before and after applying the TFPF method for an arbitrarily selected single recording. The right-hand column shows the spectra of the filtered signals. The fitted parameters of the monoexponential synthetic signals were estimated, and the SNR and MAPE were also calculated in each case, as summarized in Table 2. The estimated parameters had large fitting errors after stacking. For example, the SNRs were low, and the MRS signal parameter T2* was overestimated in all cases. When the TFPF method was used, the retrieved model parameters from the single recordings were more consistent with the synthetic model parameters. A high increase in SNR was gained, and the MAPE largely decreased, which showed that the evaluating estimation was more accurate. Table 2. Comparison of the single recording processed by TFPF and stacking. The estimated parameters, SNR and MAPE are reported based on the synthetic MRS signal with E0 = 75 nV, T2* = 100 ms, Δf = 0 Hz and φ0 = −0.53 rad.   E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  # site 1:Q1  16 stacks  73.6 ± 3.4  133.4 ± 12.7  −0.46 ± 0.41  0.31  1.22  139.59 ± 31.25  1 stack + TFPF  74.1 ± 2.3  102.6 ± 3.0  0.03 ± 0.04  −0.51  22.6  9.32 ± 3.05  # site 1:Q2  16 stacks  76.8 ± 4.7  168.6 ± 14.5  0.43 ± 0.45  0.43  −0.74  143.5 ± 34.8  1 stack + TFPF  76.4 ± 2.6  103.3 ± 4.2  −0.04 ± 0.05  −0.48  27.48  5.74 ± 3.23  # site 2:Q1  32 stacks  69.9 ± 8.9  876.0 ± 22.4  0.91 ± 0.59  −1.75  −4.7  322.29 ± 70.06  1 stack + TFPF  77.5 ± 3.7  106.8 ± 5.5  −0.08 ± 0.12  −0.62  11.3  27.2 ± 6.2  # site 2:Q2  32 stacks  62.9 ± 11.4  835.2 ± 27.3  1.25 ± 0.55  1.12  −4.51  283.88 ± 65.65  1 stack + TFPF  71.9 ± 4.1  107.1 ± 6.3  0.09 ± 0.11  −0.65  12.3  29.28 ± 6.50    E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  # site 1:Q1  16 stacks  73.6 ± 3.4  133.4 ± 12.7  −0.46 ± 0.41  0.31  1.22  139.59 ± 31.25  1 stack + TFPF  74.1 ± 2.3  102.6 ± 3.0  0.03 ± 0.04  −0.51  22.6  9.32 ± 3.05  # site 1:Q2  16 stacks  76.8 ± 4.7  168.6 ± 14.5  0.43 ± 0.45  0.43  −0.74  143.5 ± 34.8  1 stack + TFPF  76.4 ± 2.6  103.3 ± 4.2  −0.04 ± 0.05  −0.48  27.48  5.74 ± 3.23  # site 2:Q1  32 stacks  69.9 ± 8.9  876.0 ± 22.4  0.91 ± 0.59  −1.75  −4.7  322.29 ± 70.06  1 stack + TFPF  77.5 ± 3.7  106.8 ± 5.5  −0.08 ± 0.12  −0.62  11.3  27.2 ± 6.2  # site 2:Q2  32 stacks  62.9 ± 11.4  835.2 ± 27.3  1.25 ± 0.55  1.12  −4.51  283.88 ± 65.65  1 stack + TFPF  71.9 ± 4.1  107.1 ± 6.3  0.09 ± 0.11  −0.65  12.3  29.28 ± 6.50  View Large 5 FIELD EXAMPLES In this section, we showed the de-noising results, parameter estimations and MRS inversion results, which were associated with the processed signals from an actual MRS survey. The data, which could be used to compare the results of the stacking and TFPF methods, were obtained from the regional groundwater investigations near the city of Changchun, Jilin Province, China. The Earth's magnetic field had an intensity of 54 392 nT. All measurements were conducted in Shaoguo using the JLMRS-I system developed by Jilin University. The observed datasets were obtained by using a 100 m side square loop with one turn. Sixteen pulse moments of 224.4–8554.8 A ms were used. Fig. 7 shows the initial data and processed data by the stacking and TFPF methods. Because the initial data were intensely noisy as shown in Fig. 7(a), a pre-processing procedure was preliminarily used to remove the spikes and power-line harmonics, remaining the residual random noise. Then, the result was further processed by applying the traditional stacking method, as shown in Fig. 7(b). In addition, an attempt was made to suppress the random noise of the single recording, and Fig. 7(c) shows the de-noising result. Figure 7. View largeDownload slide Field example: (a) initial data consist of spikes and power-line harmonics from an actual MRS survey in the suburban areas of Changchun, Jilin, China. The processed result by averaging all stacks (b) and single recording processed by TFPF method (c). The initial data and processed results are shown as functions of the measurement time and pulse moment. Figure 7. View largeDownload slide Field example: (a) initial data consist of spikes and power-line harmonics from an actual MRS survey in the suburban areas of Changchun, Jilin, China. The processed result by averaging all stacks (b) and single recording processed by TFPF method (c). The initial data and processed results are shown as functions of the measurement time and pulse moment. Fig. 8 shows the de-noising results of the measurement recordings from four of sixteen pulse moments, which contained both signal and noise. The left-hand column was processed by the stacking method, where 16 stacks were used. The decaying exponential form was observed in Figs 8(a) and (d), but the SNRs remained low in Figs 8(g) and (j). The middle column shows 1 stacked signal after the TFPF method. This achieved better results than stacking, even for the complex noise condition in Figs 8(h) and (k). The spectral analysis in the right-hand column also showed that the TFPF method was superior to stacking. From the inset graphs in panels (i) and (l), we could see that the harmonic frequency was still present after pre-processing. Additionally, the real noise recording was not pure white, but the TFPF method performed well in these cases. Figure 8. View largeDownload slide Field example: time-series and spectra of four of sixteen pulse moments, which contained both filtered signal and noise in Changchun, Jilin, China. The left-hand column is 16 stacked signal and noise. The middle column is 1 stacked signal after TFPF method. The right-hand column shows the spectra of filtered signals. The inset graphs in panel (i) and (l) zoom in the spectra. Figure 8. View largeDownload slide Field example: time-series and spectra of four of sixteen pulse moments, which contained both filtered signal and noise in Changchun, Jilin, China. The left-hand column is 16 stacked signal and noise. The middle column is 1 stacked signal after TFPF method. The right-hand column shows the spectra of filtered signals. The inset graphs in panel (i) and (l) zoom in the spectra. Parameter estimations for these soundings against the pulse moments are shown in Fig. 9. It should be noted that the difference in initial amplitude between stacking and TFPF was not significant, but the relaxation time obtained by TFPF was more accurate than that by stacking. Because the estimated frequency offset corresponded to zero and the phase varied smoothly, it could be concluded that the magnetic resonance signal was reliably detected. Figure 9. View largeDownload slide Field example: parameter estimations of 16 stacked signals and 1 stacked signal processed by the TFPF method in Changchun, Jilin, China. Figure 9. View largeDownload slide Field example: parameter estimations of 16 stacked signals and 1 stacked signal processed by the TFPF method in Changchun, Jilin, China. The MRS inversion results, together with the corresponding borehole log, are shown in Fig. 10. The inversion results obtained using the QT inversion method (Müller-Petke & Yaramanci 2010) that corresponded to the water content and transverse relaxation time were illustrated in the upper and lower rows, respectively. From the inversion interpretation, the layer boundaries were consistent with the borehole log. Two main aquifers were detected: a fine to coarse sand stratum at shallow depth (5–20 m) and a fine and middle sand stratum at depths of 37.5–50 m. The middle mudstone stratum of no water from 39–45 m could not be detected because of the vertical resolution. From the inversion result, which was associated with the water content, the aquifer that was detected by the TFPF method at 37.5–39 m was more accurate. Good estimates of the transverse relaxation time were also obtained. However, the performance in the deep layer was not good because of the poor SNRs corrupted by random noise. Moreover, the inversion results in Fig. 10 were from a real field site with better quality data. If the noise obtained from other sites had complex noise disturbance, we agree that the TFPF method would obtain better results. From a practical perspective and for a quick measurement, the procedure of fewer stacks using the TFPF method is efficient. Figure 10. View largeDownload slide Field example: borehole log and MRS inversion results with the stacking method (16 stacks) and TFPF method (1 stack signal) in Changchun, Jilin, China. The upper row denotes the water content, and the lower row is the relaxation time. Figure 10. View largeDownload slide Field example: borehole log and MRS inversion results with the stacking method (16 stacks) and TFPF method (1 stack signal) in Changchun, Jilin, China. The upper row denotes the water content, and the lower row is the relaxation time. 7 CONCLUSION In the present study, we investigated the TFPF method for attenuating the random noise of MRS signals. Based on the theory of instantaneous frequency estimation and noise accumulation property, we found that the TFPF method effectively attenuated the noise and preserved the signal. We tested the method on simulated data and estimated the signal parameters, and the results showed that the TFPF method could work well at a noise level with an SNR as low as −17 dB. To further demonstrate the effectiveness of our method, the real MRS noise with synthetic MRS signals were used, which showed that the TFPF method could effectively mitigate random noise and preserve the signal component. Finally, the TFPF method was used to process the real MRS data and showed that the TFPF method performed better than the traditional stacking method for 16 stacks. Our results suggest that if stacking is an essential procedure for signal processing, the number of stacks can be appropriately reduced to shorten the measurement time and improve the measurement efficiency using the TFPF technology. Nevertheless, the type of dataset is the limitation of the present method. One effective way of applying the TFPF method to process the oscillating MRS signal is increasing the linearity of the signal. Demodulating the oscillating signal will help increase the linearity, but the amount of the dataset is large, and the runtime is slow when processing the demodulated signal with the TFPF. In the future, the TFPF should be combined with other methods to increase its manipulability in the oscillating data. Acknowledgements We would like to thank Mike Müller-Petke for fruitful discussions on the TFPF method. We would also like to thank the editor and reviewers for their careful and helpful comments which help us to improve our manuscript. This work was supported in part by the National Foundation of China under Grants 41722405, 41374075, 41604095 and 41504086, the National Key Foundation for Exploring Scientific Instrument under Grant 2011YQ030133, Graduate Innovation Fund of Jilin University under Project 2017003, and the Jilin Outstanding Professor Plan 20150519008JH. REFERENCES Arnold M.J., Boashash B., 1994. Time-frequency peak filtering: analysis and implementation details, in Proceedings of the IEEE-SP International Symposium on Time-Frequency & Time-Scale Analysis , Philadelphia, PA, 256– 259. Arnold M.J., Roessgen M., Boashash B., 1994. Filtering real signals through frequency modulation and peak detection in the time-frequency plane, in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’94) , Adelaide, SA, Australia, pp. 345– 348. Behroozmand A.A., Dalgaard E., Christiansen A.V., Auken E., 2013. A comprehensive study of parameter determination in a joint MRS and TEM data analysis scheme, Near Surf. Geophys. , 11: 557– 567. https://doi.org/10.3997/1873-0604.2013040 Google Scholar CrossRef Search ADS   Behroozmand A.A., Keating K., Auken E., 2015. A review of the principles and applications of the NMR technique for near-surface characterization, Surv. Geophys. , 36( 1), 27– 85. https://doi.org/10.1007/s10712-014-9304-0 Google Scholar CrossRef Search ADS   Bernard J., 2006. Instruments and field work to measure a magnetic Resonance Sounding with NUMIS equipment, in 3rd Magnetic Resonance Sounding Workshop, Madrid, Spain. Boashash B., 1992. Estimating and interpreting the instantaneous frequency of a signal—part 1: fundamentals; part 2: algorithms and applications, Proc. IEEE , 80( 4), 519– 568. Boashash B., 2003. Time-Frequency Signal Analysis and Processing , 2nd edn, pp. 663– 670, Elsevier Press. Boashash B., Mesbah M., 2004. Signal enhancement by time-frequency peak filtering, IEEE Trans. Signal Process. , 52( 4), 929– 937. https://doi.org/10.1109/TSP.2004.823510 Google Scholar CrossRef Search ADS   Dalgaard E., Auken E., Larsen J., 2012. Adaptive noise cancelling of multichannel magnetic resonance sounding signals, Geophys. J. Int. , 191( 1), 88– 100. https://doi.org/10.1111/j.1365-246X.2012.05618.x Google Scholar CrossRef Search ADS   Dlugosch R., Günther T., Lukàcs T., Müller-Petke M., 2016. Localization and identification of thin oil layers using a slim-borehole nuclear magnetic resonance tool, Geophysics , 81( 4), WB109– WB118. https://doi.org/10.1190/geo2015-0464.1 Google Scholar CrossRef Search ADS   Gao J.Z., 2011. Detection of Weak Signals , 2nd edn, Beijing, Tsinghua University Press. Goswami J.C., Hoefel A.E., 2004. Algorithms for estimating instantaneous frequency, Signal Process. , 84( 8), 1423– 1427. https://doi.org/10.1016/j.sigpro.2004.05.016 Google Scholar CrossRef Search ADS   Grunewald E., Walsh D., 2013. Multiecho scheme advances surface NMR for aquifer characterization, Geophys. Res. Lett. , 40( 24), 6346– 6350. https://doi.org/10.1002/2013GL057607 Google Scholar CrossRef Search ADS   Grunewald E., Grombacher D., Walsh D., 2016. Adiabatic pulses enhance surface nuclear magnetic resonance measurement and survey speed for groundwater investigations, Geophysics , 81( 4), WB85– WB96. https://doi.org/10.1190/geo2015-0527.1 Google Scholar CrossRef Search ADS   Hertrich M., 2008. Imaging of groundwater with nuclear magnetic resonance, Prog. Nucl. Magn. Reson. Spectrosc. , 53( 4), 227– 248. https://doi.org/10.1016/j.pnmrs.2008.01.002 Google Scholar CrossRef Search ADS   Hyndman R.J., Koehler A.B., 2006. Another look at measures of forecast accuracy, Int. J. Forecast. , 22( 4), 679– 688. https://doi.org/10.1016/j.ijforecast.2006.03.001 Google Scholar CrossRef Search ADS   Jiang C.D., Lin J., Duan Q.M., Sun S.Q., Tian B.F., 2011. Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements, Near Surf. Geophys. , 9( 5), 459– 468. Larsen J.J., 2016. Model-based subtraction of spikes from surface nuclear magnetic resonance data, Geophysics , 81( 4), WB1– WB8. https://doi.org/10.1190/geo2015-0442.1 Google Scholar CrossRef Search ADS   Larsen J.J., Dalgaard E., Auken E., 2014. Noise cancelling of MRS signals combining model-based removal of powerline harmonics and multichannel Wiener filtering, Geophys. J. Int. , 196( 2), 828– 836. https://doi.org/10.1093/gji/ggt422 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 1998. Processing of surface proton magnetic resonance signals using non-linear fitting, J. Appl. Geophys. , 39( 2), 77– 83. https://doi.org/10.1016/S0926-9851(98)00011-1 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 2002. A review of the basic principles for proton magnetic resonance sounding measurements, J. Appl. Geophys. , 50( 1–2), 3– 19. https://doi.org/10.1016/S0926-9851(02)00127-1 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 2003. Removal of power-line harmonics from proton magnetic resonance measurements, J. Appl. Geophys. , 53( 2–3), 103– 120. https://doi.org/10.1016/S0926-9851(03)00041-7 Google Scholar CrossRef Search ADS   Legchenko A., Baltassat J.-M., Beauce A., Berbard J., 2002. Nuclear magnetic resonance as a geophysical tool for hydrogeologists, J. Appl. Geophys. , 50( 1–2), 21– 46. https://doi.org/10.1016/S0926-9851(02)00128-3 Google Scholar CrossRef Search ADS   Müller-Petke M., Costabel S., 2014. Comparison and optimal parameter setting of reference-based harmonic noise cancellation in time and frequency domain for surface-NMR, Near Surf. Geophys. , 12( 2), 199– 210. Müller-Petke M., Yaramanci U., 2010. QT inversion – Comprehensive use of the complete surface NMR data set, Geophysics , 75( 4), WA199– WA209. https://doi.org/10.1190/1.3471523 Google Scholar CrossRef Search ADS   Roessgen M., Boashash B., 1994. Time-frequency peak filtering applied to FSK signals, in Proceedings of the IEEE-SP International Symposium on Time-Frequency & Time-Scale Analysis , Philadelphia, PA, pp. 516– 519. Strehl S., 2006. Development of Strategies for Improved Filtering and Fitting of SNMR Signals, MSc thesis, Technical University of Berlin. Walsh D. et al.  , 2013. A small-diameter NMR logging tool for groundwater investigations, Groundwater , 51( 6), 914– 926. https://doi.org/10.1111/gwat.12024 Google Scholar CrossRef Search ADS   Walsh D.O., 2008. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigations, J. Appl. Geophys. , 66( 3–4), 140– 150. https://doi.org/10.1016/j.jappgeo.2008.03.006 Google Scholar CrossRef Search ADS   Wang Z.X., 2009. The key technology of magnetic resonance sounding instrument for groundwater investigation, DSc thesis , Jilin University. Wu N., Li Y., Tian Y.N., Zhong T., 2016. Trace-transform-based time-frequency filtering for seismic signal enhancement in Northeast China, C. R. Geosci. , 348( 5), 360– 367. https://doi.org/10.1016/j.crte.2016.02.001 Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Time–frequency peak filtering for random noise attenuation of magnetic resonance sounding signal

, Volume 213 (2) – May 1, 2018
12 pages

/lp/ou_press/time-frequency-peak-filtering-for-random-noise-attenuation-of-magnetic-xM2qE7n7v0
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy001
Publisher site
See Article on Publisher Site

### Abstract

Summary When measuring in a geomagnetic field, the method of magnetic resonance sounding (MRS) is often limited because of the notably low signal-to-noise ratio (SNR). Most current studies focus on discarding spiky noise and power-line harmonic noise cancellation. However, the effects of random noise should not be underestimated. The common method for random noise attenuation is stacking, but collecting multiple recordings merely to suppress random noise is time-consuming. Moreover, stacking is insufficient to suppress high-level random noise. Here, we propose the use of time–frequency peak filtering for random noise attenuation, which is performed after the traditional de-spiking and power-line harmonic removal method. By encoding the noisy signal with frequency modulation and estimating the instantaneous frequency using the peak of the time–frequency representation of the encoded signal, the desired MRS signal can be acquired from only one stack. The performance of the proposed method is tested on synthetic envelope signals and field data from different surveys. Good estimations of the signal parameters are obtained at different SNRs. Moreover, an attempt to use the proposed method to handle a single recording provides better results compared to 16 stacks. Our results suggest that the number of stacks can be appropriately reduced to shorten the measurement time and improve the measurement efficiency. Hydrogeophysics, Fourier analysis, Time-series analysis 1 INTRODUCTION Magnetic resonance sounding (MRS) is a geophysical method that can directly detect water because of its unique sensitivity to protons, which can be used to provide the depth profile of water content in the subsurface (e.g. Legchenko & Valla 2002; Hertrich 2008; Behroozmand et al.2015). MRS has been widely applied in hydrogeological investigations (Legchenko et al.2002). In recent years, this method has made rapid developments in logging (e.g. Walsh et al.2013), laboratory studies (e.g. Dlugosch et al.2016), and field measurements (e.g. Grunewald et al.2016). However, when detecting in the geomagnetic field (∼0.05 mT), the weak MRS signals are often corrupted by noise, which hinders the estimation of important parameters and restricts a more widespread use of the MRS method. The noise interferences that disturb the MRS signals mainly include spiky noise, power-line harmonic noise and random noise. The common de-noising workflow is as follows: (1) de-spiking, (2) harmonic noise cancellation and (3) random noise attenuation. To improve the signal-to-noise ratio (SNR) and better estimate the signal parameters, several effective noise-cancelling approaches are proposed, which are focused mainly on de-spiking and harmonic noise cancellation. Jiang et al. (2011) suggested a statistical stacking method to discard spikes based on the Romanovsky criterion. Dalgaard et al. (2012) implemented a spike detection algorithm based on a nonlinear energy operator. A spike from an electric fence was modelled in a study by Larsen (2016) as an impulsive excitation of a tuneable fourth-order bandpass filter and subtracted. Notch filtering (Legchenko & Valla 2003) was the most common and effective method for removing power-line harmonics in single-channel instruments. Then, based on the multichannel instrument presented by Walsh (2008), advanced techniques for noise cancellation were enabled and several methods were proposed (e.g. Walsh 2008; Dalgaard et al.2012; Müller-Petke & Costabel 2014). When there were multiple noise sources, a de-noising strategy that cancelled the power-line harmonics through a model-based approach combined with multichannel Wiener filtering was suggested by Larsen et al. (2014). It should be noted that both spiky noise from electric fences and power-line harmonic noise have specific mathematical expressions and can be methodically removed. Random noise commonly includes white noise, coloured noise, band-limited white noise, narrowband noise, etc. White noise is most common in electronic circuits with a constant power spectral density and homogeneous frequency components. In contrast, coloured noise with a wide band has a non-uniform power spectral density function (Gao 2011). Band-limited white noise and narrowband noise can be considered as white noise convolved with the impulse responses of an ideal low-pass filter and an ideal bandpass filter, respectively. The power spectral density of band-limited white noise is constant within a limited bandwidth of low frequency, and the power spectral density function of narrowband noise is limited to a narrow bandwidth. Random noise cannot be described by a specific mathematical expression, as it is difficult to remove. The methods involving random noise attenuation have rarely been reported in the current literature. At present, the most widely used method is stacking, that is, averages of multiple measurements (e.g. Legchenko & Valla 2002; Dalgaard et al.2012). However, the effect of stacking is limited in a noisy environment, and stacking is time-consuming. In the present study, we propose the use of time–frequency peak filtering (TFPF; Boashash 2003) for random noise attenuation in MRS. TFPF is a new method based on the theory of instantaneous frequency estimation and the property of noise accumulation. It was first developed by Arnold and Roessgen (Arnold & Boashash 1994; Arnold et al.1994; Roessgen & Boashash 1994) and was applied to newborn electroencephalograms (e.g. Boashash & Mesbah 2004) and seismic fields (e.g. Wu et al.2016). In the theoretical derivation part of this paper, the noisy signal is first transformed into the instantaneous frequency (IF) of an analytic signal, and the analysis of the analytic signal will then produce a time–frequency distribution and the underlying MRS signal is finally achieved by taking the peak of the time–frequency distribution. The theoretical considerations were validated using numerical simulations, and the signal parameters were well estimated at different SNRs. In addition, the method was tested on a synthetic MRS signal contaminated with real noise recordings, and the noise recordings were collected from two different field surveys. The test showed that the attempt to use the proposed method to handle a single recording provided a better result than stacking. Furthermore, the real field testing results were provided and discussed. We showed that random noise could be effectively suppressed, and an improved retrieval of the underlying MRS signal was acquired. 2 RANDOM NOISE IN MRS RECORDINGS The signals to be processed in this paper are the envelopes provided by MRS devices. The basic introduction to the envelope recordings is adequately described by Legchenko & Valla (2002). The original MRS signal is generated by the excited protons oscillating with the Larmor frequency, which is an exponentially decaying function of time as follows:   $$E(t,q) = {E_0}(q) \cdot {e^{{{ - t} / {T_2^*(q)}}}} \cdot \cos (2\pi {f_L}t + {\varphi _0}(q)),$$ (1)where t is the time sample, E0 is the initial amplitude of the MRS signal and is proportional to the water content, $$T_2^*$$ is the transverse relaxation time correlated to the pore geometry in an aquifer, fL is the Larmor frequency and may slightly vary from one place to another, φ0 is the phase shift between the recorded signal and the exciting current, q is the energizing pulse moment. The actual MRS signal is a multiexponential decay, but eq. (1) only explains a monoexponential, which can be considered an oversimplification. A high sampling frequency is required to measure the oscillating MRS signal, that is, dozens of kHz, so the amount of recorded data is large. To reduce the requirement for data acquisition and increase the measurement efficiency, a two-channel MRS signal detector followed by a low-pass filter is usually used to measure the envelope signal and the MRS devices, which record the envelope signal include JLMRS-I system (Jiang et al.2011), NumisPlus and NumisLite (Bernard 2006),   \begin{eqnarray} {E_{{\rm{Detection}}}}(t) &=& {E_0} \cdot {e^{{{ - t} / {T_2^*}}}} \cdot \cos (2\pi \Delta ft + {\varphi _0}) + {\varepsilon _X}(t)\nonumber\\ &&+\, i \cdot \left( {{E_0} \cdot {e^{{{ - t} / {T_2^*}}}} \cdot \sin (2\pi \Delta ft + {\varphi _0}) + {\varepsilon _Y}(t)} \right), \end{eqnarray} (2)where EDetection(t) denotes the envelope signal and Δf is the frequency offset of the Larmor frequency and the frequency of the excitation field. It is desirable that the frequency offset approach zero as much as possible. εX(t) and εY(t) are the additive measurement noises superimposed on the in phase and out phase of the MRS signal, which mainly consist of spiky noise, power-line harmonic noise and random noise. The random noise in the MRS originates from wide band background noise and receiver electronics noise. A field case of random noise attenuation by stacking is illustrated in Fig. 1. Figure 1. View largeDownload slide Field case: time-series and spectra of 1, 4, 8 and 16 stacks at one pulse moment, respectively. The left-hand column is the time-series. The blue lines display signals during the MRS experiments. The dashed black lines display noise recorded prior to the MRS signal measurement. The right-hand column is the spectra. The blue lines display the spectra of the signals, and the dashed black lines display the spectra of the noise that correspond to the left-hand column. Figure 1. View largeDownload slide Field case: time-series and spectra of 1, 4, 8 and 16 stacks at one pulse moment, respectively. The left-hand column is the time-series. The blue lines display signals during the MRS experiments. The dashed black lines display noise recorded prior to the MRS signal measurement. The right-hand column is the spectra. The blue lines display the spectra of the signals, and the dashed black lines display the spectra of the noise that correspond to the left-hand column. Fig. 1 shows the measurement results from 1, 4, 8, and 16 stacks at one pulse moment, which contained both signals and noise. The Larmor frequency was 2326 Hz. As widely used MRS working schemes, the field noise was directly recorded prior to the signal measurement. Here, spikes and power-line harmonics were removed prior to stacking. The time-series of the signal was corrupted by noise, and the frequency response of the signal was almost invisible from the spectral analysis for 1 stacked signal, as shown in Figs 1(a) and (b). Figs 1(c) and (d) illustrate that the noise was suppressed, and the signal was enhanced when the number of stacks was 4. When the number of stacks continued to increase, the noise level decreased. The comparisons between 8 stacks (row 3) and 16 stacks (row 4) showed that the SNR increased with the increasing stacks, but the random noise still remained. Assuming that si is the signal of the single recording and εi is the noise of the single recording, the SNR after N stacks is given as follows:   $${\rm {SNR}} = \frac{{\sum\limits_{i = 1}^{{N_s}} {{s_i}} }}{{\sqrt {\sum\limits_{i = 1}^{{N_s}} {{\varepsilon _i}^2} } }} = \frac{{{N_s} \cdot s}}{{\sqrt {{N_s}} \cdot \varepsilon }} = \sqrt {{N_s}} \cdot \frac{s}{\varepsilon },$$ (3)where $$\sum\limits_{i = 1}^{{N_s}} {{s_i}}$$ and $$\sqrt {\sum\limits_{i = 1}^{{N_s}} {{\varepsilon _i}^2} }$$ are the sum of the signals and noise intensity for Ns stacks, respectively; s and ε are the mean magnitude of the signals and noise for Ns stacks, respectively; $$\frac{s}{\varepsilon }$$ is the SNR before stacking, so that the SNR increases $$\sqrt {{N_s}}$$ times (where Ns is the number of multiple measurements). This relationship is obtained and shown in Fig. 2. Figure 2. View largeDownload slide Schematic of the signal-to-noise ratio with increasing number of stacks. Figure 2. View largeDownload slide Schematic of the signal-to-noise ratio with increasing number of stacks. Fig. 2 shows the schematic of the SNR with the increased number of stacks. This observation implies that the SNR cannot always be increased rapidly with the number of stacks and it has diminishing returns. From 1 stack to approximately 32 stacks, the SNR notably and rapidly increases; from 32 stacks to 100 stacks, it steadily increases. Then, the SNR increases very slowly after 100 stacks; that is, increasing the stack numbers is meaningless because the measurement takes a very long time. Here, we take 100 stacks as an example. If the multipulse Carr–Purcell–Meiboom–Gill (CPMG) sequence is used in the MRS measurement (Grunewald & Walsh 2013) and a full set of phase cycles is composed of four individual recordings, a complete CPMG data set takes 20 s. To obtain the accurate hydrologic information, 16 pulse moments need to be carried out. If each full cycle were repeated 100 times for noise reduction, the overall measurement time would be 12 or 13 h, including the time required to move the equipment and place the measurement coils in the field. Therefore, collecting multiple recordings to obtain the best SNR is time-consuming, and the method to average the time-series is not sufficient when the proton magnetic resonance response is covered by random noise. A more effective method is required. Thus, we developed the TFPF method, which handles a single recording to replace the stacking method for random noise attenuation. 3 METHODOLOGY 3.1 Time–frequency peak filtering algorithm Implementation of the TFPF algorithm is simple, whereby the noisy signal to be filtered is first encoded as the IF of a frequency-modulated (FM) analytic signal, and the analytic signal plays a key role in the TFPF algorithm. Analysing the analytic signal with the pseudo-Wigner-Ville distribution (PWVD) will produce a time–frequency distribution exhibiting significant energy concentration around the IF of the analytic signal. The IF is then estimated by taking the peak of the time–frequency distribution to recover the desired signal without losing useful information. When implementing the TFPF algorithm with PWVD, two special steps are introduced to improve the performance, including endpoint extension and scaling. The practical application of the TFPF algorithm is illustrated in Fig. 3, and a detailed description of the TFPF scheme is presented in the following steps. Figure 3. View largeDownload slide Flowchart of the TFPF algorithm. The TFPF algorithm consists of two stages. In the first stage, the noisy signal is encoded as the instantaneous frequency of the analytic signal after applying endpoint extension and scaling. In the second stage, the time–frequency distribution of the analytic signal is calculated using the PWVD, then take the peak of the PWVD of the analytic signal, followed by IF conversion, inverse scaling and endpoint rejecting, the filtered signal without random noise interference is finally obtained. Figure 3. View largeDownload slide Flowchart of the TFPF algorithm. The TFPF algorithm consists of two stages. In the first stage, the noisy signal is encoded as the instantaneous frequency of the analytic signal after applying endpoint extension and scaling. In the second stage, the time–frequency distribution of the analytic signal is calculated using the PWVD, then take the peak of the PWVD of the analytic signal, followed by IF conversion, inverse scaling and endpoint rejecting, the filtered signal without random noise interference is finally obtained. Let x(m) be the random-noise-corrupted MRS signal, which is modelled by the following equation:   $$x(m) = s(m) + \varepsilon (m),$$ (4)where m is the discrete time sample, m = 1, 2, …, M, M is the length of x(m), s(m) is the desired MRS signal and ε(m) is the random noise. 3.1.1 First stage: encoding First step: extend the endpoints of the sampled signal. To obtain the unbiased estimation of the desired signal, an appropriate window length is selected during PWVD. However, the available sampling points at the boundaries are less than the window length, and the energy would not concentrate at the boundaries of the analytic signal's IF on the time–frequency plane. As a result, the values of the estimated signal at both ends are deviated from the real values when the observed dataset is directly processed. To avoid endpoint distortion, endpoint extension should be performed, and the extended samples are required to be larger than half of the window length. Here, we chose a simple way to extend the observed dataset with a length of p samples at both ends, which is given as follows:   $$x'(n) = \left\{ \begin{array}{@{}l@{}} x(1),1 \le n \le p\\ x(n - p),p + 1 \le n \le p + M\\ x(M),p + M + 1 \le n \le 2p + M \end{array} \right.,$$ (5)where x'(n) is the signal after endpoint extension and 1 ≤ n ≤ 2p + M. Assuming that the window length W = 2L + 1, p ≥ L + 1 is required. Second step: scale the amplitude of the extended signal. Modulating the unscaled signal will lead to aliasing, which produces discontinuities in the estimated IF at the frequency boundaries of the time–frequency plane. To process all signals with different amplitudes, scaling should be performed before encoding to provide suitable frequency limits on the encoded signal. For easy understanding, the signal which is then encoded as the IF of the analytic signal is sampled at a normalized frequency of 1 Hz, the amplitude of the signal should range from half of the normalized frequency (0.5) to 0:   $${x_c}(n) = (a - b)\frac{{x'(n) - \min [x'(n)]}}{{\max [x'(n)] - \min [x'(n)]}} + b,$$ (6)where xc(n) is the scaled signal. After applying amplitude scaling, parameters a and b are the maximum and minimum of the discrete noisy signal, respectively, which satisfies the constraint 0 ≤ b < a ≤ 0.5. To preserve the signal details to make the amplitude of the scaled signal vary as large as possible, we normally let a = 0.5 and b = 0. Third step: encode the scaled signal as the IF of an analytic signal. According to the notion of IF (e.g. Goswami & Hoefel 2004), the IF of a given analytic signal z(t) = ej2πϕ(t) is as follows:   $${f_z}(t) = \frac{1}{{2\pi }}\frac{{{\rm{d}}\phi (t)}}{{{\rm{d}}t}},$$ (7)where ϕ(t) is the instantaneous phase and fz(t) is the IF of the analytic signal, and thus, the analytic signal can also be expressed as follows:   $$z(t) = {{\rm{e}}^{j2\pi \int_{{ - \infty }}^{t}{{{f_z}(\lambda ){\rm{d}}\lambda }}}}.$$ (8) Likewise, the discrete scaled signal xc(n) is encoded as the IF of the unit-amplitude FM analytic signal z(n), given as follows:   $$z(n) = {{\rm{e}}^{j2\pi \sum\limits_{\lambda = 0}^n {{x_c}(\lambda )} }}.$$ (9) 3.1.2 Second stage: IF estimation First step: complete the PWVD of the analytic signal. Calculate the time–frequency distribution of the analytic signal using PWVD, the discrete PWVD of the analytic signal z(n) can be expressed as follows:   $${\rm{PWV}}{{\rm{D}}_z}(n,k) = 2\sum\limits_{l = - L}^L {{\rm w}(l)z(n + l){z^*}(n - l){{\rm{e}}^{ - j4\pi kl}}} ,$$ (10)where w(l) is a window function, the symbol * denotes the complex conjugate, k is the discrete frequency sample, and 2L + 1 is the window length (W) considered in the PWVD discrete time implementation. It is important that the signal be linear to obtain an unbiased estimation (Boashash & Mesbah 2004). However, the MRS signal is a nonlinear function of time, so an appropriate window length should be selected to make the signal as close to linear as possible across the window length. The relationship related to the TFPF window length is given for the case of signal estimation, which is expressed as follows:   $$W \le \frac{{0.634{f_s}}}{{\pi f}}.$$ (11) Eq.(11) gives the maximum window length as a function of sampling rate (fs) and signal frequency (f). For the signals generated by protons, the frequency of the envelope signal would not exceed the range of 5 to −5 Hz (Strehl 2006). To reduce the bias, we select a smaller window length, and the signal frequency of 5 Hz is always chosen to calculate the window length. In addition, the sampling rate is normally four times the Larmor frequency. Without a loss of signal characteristics, resampling of the time-series is always needed, for example, the resampling rate of the envelope signal is a quarter of the Larmor frequency (Wang 2009). The random noise that corrupted the envelope signals can be attenuated effectively by the TFPF method, but it cannot process the oscillating signals well at present because the oscillating NMR signals recorded at dozens of kHz sampling rates have a high nonlinearity and the IF estimations are biased. Signal preservation and noise suppression are difficult to weigh. If the selected window length is large, the signal component will be lost. Conversely, the noise will not be better suppressed. This is the limitation of the present method. Second step: take the frequency samples of the IF estimation. Take the peak of PWVDz(n, k) to estimate the IF as follows:   $${\hat{k}_z}(n) = {\arg _k}\max [{\rm{PWV}}{{\rm{D}}_z}(n,k)],$$ (12)where $${\hat{k}_z}(n)$$ is the frequency sample of the IF estimation, which can be considered a recovery of the desired signal based on the standard time–frequency peak-detection algorithm (Boashash 1992). The subscript k for $${\arg _k}$$ means a discrete frequency sample; when maximizing the PWVD over the frequency, the peak of PWVD is found for each value of k. Third step: convert the frequency samples to the IF estimation. Because the values of $${\hat{k}_z}(n)$$ are composed of the frequency samples, we should convert the estimated frequency samples as follows:   $${\hat{s}_c}(n) = ({\hat{k}_z}(n) - 1)\frac{{{f_s}}}{N},$$ (13)where $${\hat{s}_c}(n)$$ is the estimated IF after conversion; $$\frac{{{f_s}}}{N}$$ is the frequency resolution; fs is the sampling rate and N is the number of frequency bins during the discrete PWVD. Fourth step: inverse scaling of the signal after IF conversion. The signal after IF conversation ranges from a to b, the order of the desired signal is recovered by an inverse scaling operation:   \begin{eqnarray} \hat{s}'(n) = \frac{{({{\hat{s}}_c}(n) - b)(\max [x'(n)] - \min [x'(n)])}}{{a - b}} + \min [x'(n)],\nonumber\\ \end{eqnarray} (14)where $$\hat{s}'(n)$$ is the desired signal with a real amplitude. Fifth step: reject the extended endpoints. The final filtered result $$\hat{s}(m)$$ shall be obtained by rejecting the extended endpoints given as follows:   $$\hat{s}(m) = \hat{s}'(m + p),1 \le m \le M.$$ (15) 3.2 A complete signal example In what follows, we show a complete signal example to explain each step of the TFPF presented in the previous section. A synthetic MRS envelope signal with the parameters E0 = 200 nV and T2* = 150 ms was generated and had been corrupted by random noise with a standard deviation of 50 nV (the mean value of the artificial random noise is 0 nV in all the following experiments), making the SNR 7.11 dB. Then, the TFPF method was performed to mitigate the random noise, and the basic procedure is illustrated in Fig. 4. Figure 4. View largeDownload slide Example of the basic procedure for the random noise attenuation. (a) The original noisy signal. (b) Endpoint extension. (c) Amplitude scaling. (d) Time–frequency implementation. (e) IF estimation. (f) IF conversion. (g) Inverse scaling. (h) Endpoints rejecting. Figure 4. View largeDownload slide Example of the basic procedure for the random noise attenuation. (a) The original noisy signal. (b) Endpoint extension. (c) Amplitude scaling. (d) Time–frequency implementation. (e) IF estimation. (f) IF conversion. (g) Inverse scaling. (h) Endpoints rejecting. The random-noise-corrupted MRS signal is shown in Fig. 4(a). The length of the observed signal was 145 (Larmor frequency fL = 2320 Hz, the resampling rate was set as a quarter of the Larmor frequency, and the recording time was 250 ms). In the first stage, we dealt with the noisy signal by extending the endpoints according to eq. (5), a length of eight samples were extended at both ends of the original noisy signal, which was shown in Fig. 4(b) (first step). Then, signal scaling was performed before encoding according to eq. (6). Assuming that the extended signal, which was regarded as the IF of the analytic signal was sampled at a normalized frequency of 1 Hz, Fig. 4(c) shows the amplitude of the scaled signal was between 0.5 and 0, where we selected the scaling parameters a = 0.5 and b = 0 (second step). Finally, the scaled signal was encoded as the IF of a unit-amplitude FM analytic signal according to eq. (9). Thus, the analytic signal was obtained, and the analytic signal was important in the TFPF scheme (third step). In the second stage, the PWVD was first performed on the analytic signal according to eq. (10). The window length (W) used in the discrete PWVD was 15. A time–frequency distribution exhibiting significant energy concentration around the IF of the analytic signal is produced, which resembled the desired MRS signal as shown in Fig. 4(d). However, the random noise was scattered over the time–frequency plane due to its stochastic properties, and this yielded an appropriate separation between signal and noise components, as well as a good estimation without the effect of random noise (first step). Subsequently, the frequency samples were obtained by maximizing the PWVD along the frequency axis according to eq. (12). The number of frequency bins (N) we selected during PWVD was 128, so the maximum value of y-axis was 128 as shown in Fig. 4(e) (second step). Then, the estimated frequency samples were multiplied by the frequency resolution (fs/N) according to eq. (13), the converted signal (Fig. 4f) was generated with the amplitude ranged from 0.5 to 0 (third step). Subsequently, inverse scaling was conducted according to eq. (14), the real amplitude of the signal was recovered (Fig. 4g, fourth step), followed by the endpoints being rejected according to eq. (15) (fifth step), the filtered signal was finally obtained with parameters E0 = 201.8 ± 2.1 nV, T2* = 149.97 ± 1.60 ms and the SNR increased to 28.02 dB, which is shown in Fig. 4(h). 4 NUMERICAL SIMULATIONS In this section, to validate the proposed method, numerical experiments were performed to model the monoexponential decaying signals in the artificial random noise and real noise recordings. To quantify the effects of the random noise attenuation, a fitting algorithm described in Legchenko & Valla (1998) was used to estimate the signal parameters. The performance of the proposed method was also evaluated based on the SNR and mean absolute percentage error (MAPE). Here, as a scale-dependent error metric (Hyndman & Koehler 2006), the MAPE was quoted to measure the quality of preserved information in the de-noised signal. The MAPE is defined as follows:   $$MAPE\,{\rm{ = }}\,\frac{1}{M}\sum\limits_{m = 1}^M {\left| {\frac{{\hat{s}(m) - s(m)}}{{s(m)}} \times 100} \right|} ,$$ (16)where M is the length of samples, s(m) is the ideal signal and $$\hat{s}(m)$$ is the de-noised signal. 4.1 Synthetic example for artificial noise A monoexponential decaying signal with the parameters E0 = 200 nV, T2* = 150 ms, Δf = 0 Hz and φ0 = 1.05 rad were generated. To make the synthetic data comparable with field conditions, the added noise was simulated as follows (Behroozmand et al.2013):   \begin{eqnarray} x(m) = s(m) + {\varepsilon _{G(0,1)}}(m) \cdot {\left[\varepsilon _{{\rm{uni}}}^2(m) + {\left(\frac{{{\varepsilon _{{\rm{bac}}}}(m)}}{{s(m)}}\right)^2} \right]^{\frac{1}{2}}} \cdot s(m),\nonumber\\ \end{eqnarray} (17)where m is the discrete time sample, x(m) and s(m) denote the noisy signal and ideal signal, respectively. εG(0, 1)(m) is the Gaussian noise with a mean value of 0 and a standard deviation of 1; εuni(m) is uniform noise, which can be considered non-specified noise contributions such as structural noise; and εbac(m) is the background noise. Fig. 5 shows the results of the TFPF method for random noise attenuation under different noise conditions. To simulate MRS signal 1, the monoexponential decaying signal was corrupted by a Gaussian noise distribution with a standard deviation of 30 nV and a uniform relative noise of 3 per cent of the data values, which were assigned to Vbac and STDuni in eq. (17). The SNR was 10.81 dB. Then, the TFPF method was employed to process the noisy signal. Panel (c) showed the PWVD of the analytic signal, the observation implied that a significant energy concentration was produced, and the signal could be estimated by maximizing the PWVD of the analytic signal along the frequency axis. The signal, regarded as the IF of the analytic signal, was sampled at a normalized frequency of 1 Hz, and the frequency axis of the time–frequency plane ranged from 1 to 0 Hz. Figs 5(a) and (b) illustrated the processing results in the time and frequency domains, respectively. The random noise component was eliminated, and the signal component was preserved. The estimated parameters, SNR and MAPE after TFPF, are reported in Table 1. We observed that with TFPF, the SNR increased to 28.32 dB, the estimated E0 was 199.3 ± 2.1 nV, T2* was 151.2 ± 1.5 ms, Δf was −0.01 ± 0.02 Hz, φ0 was 1.09 rad, and MAPE was 4.8 ± 1.4 per cent, respectively. To test the effectiveness of this method for a high noise level, the Gaussian noise was increased to a standard deviation of 150 nV, and the uniform relative noise was 5 per cent of the data values, which made the SNR of the signal 2 −3.16 dB. The middle row showed the outcome of the TFPF method, where the SNR increased to 26.27 dB. The estimated MRS signal parameters E0, T2*, Δf and φ0 were 201.2 ± 2.8 nV, 149.99 ± 2.05 ms, −0.03 ± 0.05 Hz and 1.11 rad, respectively, and the calculated MAPE was 5.7 ± 1.8 per cent. Furthermore, a serious noise disturbance was simulated as signal 3. The ideal signal was contaminated with Gaussian noise, the standard deviation was 300 nV with a uniform relative noise of 10 per cent of the data values, and the SNR was −10.66 dB. We observed that the signal was submerged in a high noise level from the time-series in Fig. 5(g). The frequency offset of 0 Hz was hardly observed and indistinctly emerged from the spectrum, as shown in Fig. 5(h). The fidelity of the MRS signal was severely degraded by random noise. However, the de-noised signal showed that the TFPF had an excellent performance because the high-level random noise did not produce an energy concentration in the time–frequency plane, as shown in Fig. 5(i). The high-level random noise was mitigated without losing the signal information, and the SNR increased to 21.73 dB; the estimated E0 was 198.9 ± 4.3 nV, T2* was 146.7 ± 5.1 ms, Δf was 0.04 ± 0.07 Hz, φ0 was 0.94 rad, and MAPE was 9.2 ± 3.3 per cent. In addition, we increased the noise level and tested this method. The results showed that the underlying signals could be clean recovered in the noise level down to an SNR of −17 dB. Figure 5. View largeDownload slide Numerical experiments of modelling the monoexponential synthetic signal added to three different artificial random noise recordings. The simulated signals with three different SNRs (signal 1 = 10.81 dB, signal 2 = −3.16 dB and signal 3 = −10.66 dB) are shown in the upper, middle and lower rows, respectively. The time-series and spectra of the envelope signals are displayed in the left-hand and middle columns, respectively. The right-hand column shows the PWVD of the FM analytic signal. Figure 5. View largeDownload slide Numerical experiments of modelling the monoexponential synthetic signal added to three different artificial random noise recordings. The simulated signals with three different SNRs (signal 1 = 10.81 dB, signal 2 = −3.16 dB and signal 3 = −10.66 dB) are shown in the upper, middle and lower rows, respectively. The time-series and spectra of the envelope signals are displayed in the left-hand and middle columns, respectively. The right-hand column shows the PWVD of the FM analytic signal. Table 1. Estimated parameters, SNR and MAPE performance of TFPF on synthetic MRS signals superimposed in different artificial random noise recordings. The synthetic MRS signal has E0 = 200 nV, T2* = 150 ms, Δf = 0 Hz and φ0 = 1.05 rad.   E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  Signal 1 (SNR = 10.81 dB)  199.3 ± 2.1  151.2 ± 1.5  −0.01 ± 0.02  1.09  28.32  4.8 ± 1.4  Signal 2 (SNR = -3.16 dB)  201.2 ± 2.8  149.99 ± 2.05  −0.03 ± 0.05  1.11  26.27  5.7 ± 1.8  Signal 3 (SNR = -10.66 dB)  198.9 ± 4.3  146.7 ± 5.1  0.04 ± 0.07  0.94  21.73  9.2 ± 3.3    E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  Signal 1 (SNR = 10.81 dB)  199.3 ± 2.1  151.2 ± 1.5  −0.01 ± 0.02  1.09  28.32  4.8 ± 1.4  Signal 2 (SNR = -3.16 dB)  201.2 ± 2.8  149.99 ± 2.05  −0.03 ± 0.05  1.11  26.27  5.7 ± 1.8  Signal 3 (SNR = -10.66 dB)  198.9 ± 4.3  146.7 ± 5.1  0.04 ± 0.07  0.94  21.73  9.2 ± 3.3  View Large 4.2 Synthetic example for real noise In this subsection, another synthetic example was shown. A monoexponential decaying signal with the parameters E0 = 75 nV, T2* = 100 ms, Δf = 0 Hz and φ0 = −0.53 rad, given by eq. (2), were added to noise-only recordings from two different sites. The noise levels at site 1 were smaller than those at site 2; 16 stacks at one pulse moment for site 1 and 32 stacks at one pulse moment for site 2 were performed. We conducted two MRS experiments at each site and collected the noise-only recordings. The actual noise interferences included spiky noise, power-line harmonic noise and random noise. The serious spiky noise was found in almost all recordings. Hence, a de-spiking algorithm and a harmonic noise cancellation algorithm (Jiang et al.2011) were first used in these original noisy recordings. After de-spiking and harmonic noise cancellation, the random noise remained. Then, we arbitrarily selected a single recording from the processed multiple recordings and applied the TFPF method to suppress the residual random noise. We compared this result with stacking. Fig. 6 shows a comparison of the traditional stacking method and TFPF method for random noise attenuation. The panels in the left-hand column illustrated that the random noise component still remained after stacking, whereas good results were obtained by the TFPF method, which handled the single recordings. The processed results in the middle column showed that the residual random noise was suppressed, and the decaying signal curve appeared close to the ideal signal. Compared with the traditional stacking results, the proposed method enabled us to recover the desired MRS signal even at a low SNR. Figure 6. View largeDownload slide Comparison of the stacking and the TFPF method for the residual random noise attenuation. The panels in the left-hand column display the de-noising results by stacking. The middle column shows the results before and after applying the TFPF method for an arbitrarily selected single recording. The right-hand column shows the spectra of the filtered signals. Figure 6. View largeDownload slide Comparison of the stacking and the TFPF method for the residual random noise attenuation. The panels in the left-hand column display the de-noising results by stacking. The middle column shows the results before and after applying the TFPF method for an arbitrarily selected single recording. The right-hand column shows the spectra of the filtered signals. The fitted parameters of the monoexponential synthetic signals were estimated, and the SNR and MAPE were also calculated in each case, as summarized in Table 2. The estimated parameters had large fitting errors after stacking. For example, the SNRs were low, and the MRS signal parameter T2* was overestimated in all cases. When the TFPF method was used, the retrieved model parameters from the single recordings were more consistent with the synthetic model parameters. A high increase in SNR was gained, and the MAPE largely decreased, which showed that the evaluating estimation was more accurate. Table 2. Comparison of the single recording processed by TFPF and stacking. The estimated parameters, SNR and MAPE are reported based on the synthetic MRS signal with E0 = 75 nV, T2* = 100 ms, Δf = 0 Hz and φ0 = −0.53 rad.   E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  # site 1:Q1  16 stacks  73.6 ± 3.4  133.4 ± 12.7  −0.46 ± 0.41  0.31  1.22  139.59 ± 31.25  1 stack + TFPF  74.1 ± 2.3  102.6 ± 3.0  0.03 ± 0.04  −0.51  22.6  9.32 ± 3.05  # site 1:Q2  16 stacks  76.8 ± 4.7  168.6 ± 14.5  0.43 ± 0.45  0.43  −0.74  143.5 ± 34.8  1 stack + TFPF  76.4 ± 2.6  103.3 ± 4.2  −0.04 ± 0.05  −0.48  27.48  5.74 ± 3.23  # site 2:Q1  32 stacks  69.9 ± 8.9  876.0 ± 22.4  0.91 ± 0.59  −1.75  −4.7  322.29 ± 70.06  1 stack + TFPF  77.5 ± 3.7  106.8 ± 5.5  −0.08 ± 0.12  −0.62  11.3  27.2 ± 6.2  # site 2:Q2  32 stacks  62.9 ± 11.4  835.2 ± 27.3  1.25 ± 0.55  1.12  −4.51  283.88 ± 65.65  1 stack + TFPF  71.9 ± 4.1  107.1 ± 6.3  0.09 ± 0.11  −0.65  12.3  29.28 ± 6.50    E0 (nV)  T2* (ms)  Δf (Hz)  φ0(rad)  SNR (dB)  MAPE (%)  # site 1:Q1  16 stacks  73.6 ± 3.4  133.4 ± 12.7  −0.46 ± 0.41  0.31  1.22  139.59 ± 31.25  1 stack + TFPF  74.1 ± 2.3  102.6 ± 3.0  0.03 ± 0.04  −0.51  22.6  9.32 ± 3.05  # site 1:Q2  16 stacks  76.8 ± 4.7  168.6 ± 14.5  0.43 ± 0.45  0.43  −0.74  143.5 ± 34.8  1 stack + TFPF  76.4 ± 2.6  103.3 ± 4.2  −0.04 ± 0.05  −0.48  27.48  5.74 ± 3.23  # site 2:Q1  32 stacks  69.9 ± 8.9  876.0 ± 22.4  0.91 ± 0.59  −1.75  −4.7  322.29 ± 70.06  1 stack + TFPF  77.5 ± 3.7  106.8 ± 5.5  −0.08 ± 0.12  −0.62  11.3  27.2 ± 6.2  # site 2:Q2  32 stacks  62.9 ± 11.4  835.2 ± 27.3  1.25 ± 0.55  1.12  −4.51  283.88 ± 65.65  1 stack + TFPF  71.9 ± 4.1  107.1 ± 6.3  0.09 ± 0.11  −0.65  12.3  29.28 ± 6.50  View Large 5 FIELD EXAMPLES In this section, we showed the de-noising results, parameter estimations and MRS inversion results, which were associated with the processed signals from an actual MRS survey. The data, which could be used to compare the results of the stacking and TFPF methods, were obtained from the regional groundwater investigations near the city of Changchun, Jilin Province, China. The Earth's magnetic field had an intensity of 54 392 nT. All measurements were conducted in Shaoguo using the JLMRS-I system developed by Jilin University. The observed datasets were obtained by using a 100 m side square loop with one turn. Sixteen pulse moments of 224.4–8554.8 A ms were used. Fig. 7 shows the initial data and processed data by the stacking and TFPF methods. Because the initial data were intensely noisy as shown in Fig. 7(a), a pre-processing procedure was preliminarily used to remove the spikes and power-line harmonics, remaining the residual random noise. Then, the result was further processed by applying the traditional stacking method, as shown in Fig. 7(b). In addition, an attempt was made to suppress the random noise of the single recording, and Fig. 7(c) shows the de-noising result. Figure 7. View largeDownload slide Field example: (a) initial data consist of spikes and power-line harmonics from an actual MRS survey in the suburban areas of Changchun, Jilin, China. The processed result by averaging all stacks (b) and single recording processed by TFPF method (c). The initial data and processed results are shown as functions of the measurement time and pulse moment. Figure 7. View largeDownload slide Field example: (a) initial data consist of spikes and power-line harmonics from an actual MRS survey in the suburban areas of Changchun, Jilin, China. The processed result by averaging all stacks (b) and single recording processed by TFPF method (c). The initial data and processed results are shown as functions of the measurement time and pulse moment. Fig. 8 shows the de-noising results of the measurement recordings from four of sixteen pulse moments, which contained both signal and noise. The left-hand column was processed by the stacking method, where 16 stacks were used. The decaying exponential form was observed in Figs 8(a) and (d), but the SNRs remained low in Figs 8(g) and (j). The middle column shows 1 stacked signal after the TFPF method. This achieved better results than stacking, even for the complex noise condition in Figs 8(h) and (k). The spectral analysis in the right-hand column also showed that the TFPF method was superior to stacking. From the inset graphs in panels (i) and (l), we could see that the harmonic frequency was still present after pre-processing. Additionally, the real noise recording was not pure white, but the TFPF method performed well in these cases. Figure 8. View largeDownload slide Field example: time-series and spectra of four of sixteen pulse moments, which contained both filtered signal and noise in Changchun, Jilin, China. The left-hand column is 16 stacked signal and noise. The middle column is 1 stacked signal after TFPF method. The right-hand column shows the spectra of filtered signals. The inset graphs in panel (i) and (l) zoom in the spectra. Figure 8. View largeDownload slide Field example: time-series and spectra of four of sixteen pulse moments, which contained both filtered signal and noise in Changchun, Jilin, China. The left-hand column is 16 stacked signal and noise. The middle column is 1 stacked signal after TFPF method. The right-hand column shows the spectra of filtered signals. The inset graphs in panel (i) and (l) zoom in the spectra. Parameter estimations for these soundings against the pulse moments are shown in Fig. 9. It should be noted that the difference in initial amplitude between stacking and TFPF was not significant, but the relaxation time obtained by TFPF was more accurate than that by stacking. Because the estimated frequency offset corresponded to zero and the phase varied smoothly, it could be concluded that the magnetic resonance signal was reliably detected. Figure 9. View largeDownload slide Field example: parameter estimations of 16 stacked signals and 1 stacked signal processed by the TFPF method in Changchun, Jilin, China. Figure 9. View largeDownload slide Field example: parameter estimations of 16 stacked signals and 1 stacked signal processed by the TFPF method in Changchun, Jilin, China. The MRS inversion results, together with the corresponding borehole log, are shown in Fig. 10. The inversion results obtained using the QT inversion method (Müller-Petke & Yaramanci 2010) that corresponded to the water content and transverse relaxation time were illustrated in the upper and lower rows, respectively. From the inversion interpretation, the layer boundaries were consistent with the borehole log. Two main aquifers were detected: a fine to coarse sand stratum at shallow depth (5–20 m) and a fine and middle sand stratum at depths of 37.5–50 m. The middle mudstone stratum of no water from 39–45 m could not be detected because of the vertical resolution. From the inversion result, which was associated with the water content, the aquifer that was detected by the TFPF method at 37.5–39 m was more accurate. Good estimates of the transverse relaxation time were also obtained. However, the performance in the deep layer was not good because of the poor SNRs corrupted by random noise. Moreover, the inversion results in Fig. 10 were from a real field site with better quality data. If the noise obtained from other sites had complex noise disturbance, we agree that the TFPF method would obtain better results. From a practical perspective and for a quick measurement, the procedure of fewer stacks using the TFPF method is efficient. Figure 10. View largeDownload slide Field example: borehole log and MRS inversion results with the stacking method (16 stacks) and TFPF method (1 stack signal) in Changchun, Jilin, China. The upper row denotes the water content, and the lower row is the relaxation time. Figure 10. View largeDownload slide Field example: borehole log and MRS inversion results with the stacking method (16 stacks) and TFPF method (1 stack signal) in Changchun, Jilin, China. The upper row denotes the water content, and the lower row is the relaxation time. 7 CONCLUSION In the present study, we investigated the TFPF method for attenuating the random noise of MRS signals. Based on the theory of instantaneous frequency estimation and noise accumulation property, we found that the TFPF method effectively attenuated the noise and preserved the signal. We tested the method on simulated data and estimated the signal parameters, and the results showed that the TFPF method could work well at a noise level with an SNR as low as −17 dB. To further demonstrate the effectiveness of our method, the real MRS noise with synthetic MRS signals were used, which showed that the TFPF method could effectively mitigate random noise and preserve the signal component. Finally, the TFPF method was used to process the real MRS data and showed that the TFPF method performed better than the traditional stacking method for 16 stacks. Our results suggest that if stacking is an essential procedure for signal processing, the number of stacks can be appropriately reduced to shorten the measurement time and improve the measurement efficiency using the TFPF technology. Nevertheless, the type of dataset is the limitation of the present method. One effective way of applying the TFPF method to process the oscillating MRS signal is increasing the linearity of the signal. Demodulating the oscillating signal will help increase the linearity, but the amount of the dataset is large, and the runtime is slow when processing the demodulated signal with the TFPF. In the future, the TFPF should be combined with other methods to increase its manipulability in the oscillating data. Acknowledgements We would like to thank Mike Müller-Petke for fruitful discussions on the TFPF method. We would also like to thank the editor and reviewers for their careful and helpful comments which help us to improve our manuscript. This work was supported in part by the National Foundation of China under Grants 41722405, 41374075, 41604095 and 41504086, the National Key Foundation for Exploring Scientific Instrument under Grant 2011YQ030133, Graduate Innovation Fund of Jilin University under Project 2017003, and the Jilin Outstanding Professor Plan 20150519008JH. REFERENCES Arnold M.J., Boashash B., 1994. Time-frequency peak filtering: analysis and implementation details, in Proceedings of the IEEE-SP International Symposium on Time-Frequency & Time-Scale Analysis , Philadelphia, PA, 256– 259. Arnold M.J., Roessgen M., Boashash B., 1994. Filtering real signals through frequency modulation and peak detection in the time-frequency plane, in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’94) , Adelaide, SA, Australia, pp. 345– 348. Behroozmand A.A., Dalgaard E., Christiansen A.V., Auken E., 2013. A comprehensive study of parameter determination in a joint MRS and TEM data analysis scheme, Near Surf. Geophys. , 11: 557– 567. https://doi.org/10.3997/1873-0604.2013040 Google Scholar CrossRef Search ADS   Behroozmand A.A., Keating K., Auken E., 2015. A review of the principles and applications of the NMR technique for near-surface characterization, Surv. Geophys. , 36( 1), 27– 85. https://doi.org/10.1007/s10712-014-9304-0 Google Scholar CrossRef Search ADS   Bernard J., 2006. Instruments and field work to measure a magnetic Resonance Sounding with NUMIS equipment, in 3rd Magnetic Resonance Sounding Workshop, Madrid, Spain. Boashash B., 1992. Estimating and interpreting the instantaneous frequency of a signal—part 1: fundamentals; part 2: algorithms and applications, Proc. IEEE , 80( 4), 519– 568. Boashash B., 2003. Time-Frequency Signal Analysis and Processing , 2nd edn, pp. 663– 670, Elsevier Press. Boashash B., Mesbah M., 2004. Signal enhancement by time-frequency peak filtering, IEEE Trans. Signal Process. , 52( 4), 929– 937. https://doi.org/10.1109/TSP.2004.823510 Google Scholar CrossRef Search ADS   Dalgaard E., Auken E., Larsen J., 2012. Adaptive noise cancelling of multichannel magnetic resonance sounding signals, Geophys. J. Int. , 191( 1), 88– 100. https://doi.org/10.1111/j.1365-246X.2012.05618.x Google Scholar CrossRef Search ADS   Dlugosch R., Günther T., Lukàcs T., Müller-Petke M., 2016. Localization and identification of thin oil layers using a slim-borehole nuclear magnetic resonance tool, Geophysics , 81( 4), WB109– WB118. https://doi.org/10.1190/geo2015-0464.1 Google Scholar CrossRef Search ADS   Gao J.Z., 2011. Detection of Weak Signals , 2nd edn, Beijing, Tsinghua University Press. Goswami J.C., Hoefel A.E., 2004. Algorithms for estimating instantaneous frequency, Signal Process. , 84( 8), 1423– 1427. https://doi.org/10.1016/j.sigpro.2004.05.016 Google Scholar CrossRef Search ADS   Grunewald E., Walsh D., 2013. Multiecho scheme advances surface NMR for aquifer characterization, Geophys. Res. Lett. , 40( 24), 6346– 6350. https://doi.org/10.1002/2013GL057607 Google Scholar CrossRef Search ADS   Grunewald E., Grombacher D., Walsh D., 2016. Adiabatic pulses enhance surface nuclear magnetic resonance measurement and survey speed for groundwater investigations, Geophysics , 81( 4), WB85– WB96. https://doi.org/10.1190/geo2015-0527.1 Google Scholar CrossRef Search ADS   Hertrich M., 2008. Imaging of groundwater with nuclear magnetic resonance, Prog. Nucl. Magn. Reson. Spectrosc. , 53( 4), 227– 248. https://doi.org/10.1016/j.pnmrs.2008.01.002 Google Scholar CrossRef Search ADS   Hyndman R.J., Koehler A.B., 2006. Another look at measures of forecast accuracy, Int. J. Forecast. , 22( 4), 679– 688. https://doi.org/10.1016/j.ijforecast.2006.03.001 Google Scholar CrossRef Search ADS   Jiang C.D., Lin J., Duan Q.M., Sun S.Q., Tian B.F., 2011. Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements, Near Surf. Geophys. , 9( 5), 459– 468. Larsen J.J., 2016. Model-based subtraction of spikes from surface nuclear magnetic resonance data, Geophysics , 81( 4), WB1– WB8. https://doi.org/10.1190/geo2015-0442.1 Google Scholar CrossRef Search ADS   Larsen J.J., Dalgaard E., Auken E., 2014. Noise cancelling of MRS signals combining model-based removal of powerline harmonics and multichannel Wiener filtering, Geophys. J. Int. , 196( 2), 828– 836. https://doi.org/10.1093/gji/ggt422 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 1998. Processing of surface proton magnetic resonance signals using non-linear fitting, J. Appl. Geophys. , 39( 2), 77– 83. https://doi.org/10.1016/S0926-9851(98)00011-1 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 2002. A review of the basic principles for proton magnetic resonance sounding measurements, J. Appl. Geophys. , 50( 1–2), 3– 19. https://doi.org/10.1016/S0926-9851(02)00127-1 Google Scholar CrossRef Search ADS   Legchenko A., Valla P., 2003. Removal of power-line harmonics from proton magnetic resonance measurements, J. Appl. Geophys. , 53( 2–3), 103– 120. https://doi.org/10.1016/S0926-9851(03)00041-7 Google Scholar CrossRef Search ADS   Legchenko A., Baltassat J.-M., Beauce A., Berbard J., 2002. Nuclear magnetic resonance as a geophysical tool for hydrogeologists, J. Appl. Geophys. , 50( 1–2), 21– 46. https://doi.org/10.1016/S0926-9851(02)00128-3 Google Scholar CrossRef Search ADS   Müller-Petke M., Costabel S., 2014. Comparison and optimal parameter setting of reference-based harmonic noise cancellation in time and frequency domain for surface-NMR, Near Surf. Geophys. , 12( 2), 199– 210. Müller-Petke M., Yaramanci U., 2010. QT inversion – Comprehensive use of the complete surface NMR data set, Geophysics , 75( 4), WA199– WA209. https://doi.org/10.1190/1.3471523 Google Scholar CrossRef Search ADS   Roessgen M., Boashash B., 1994. Time-frequency peak filtering applied to FSK signals, in Proceedings of the IEEE-SP International Symposium on Time-Frequency & Time-Scale Analysis , Philadelphia, PA, pp. 516– 519. Strehl S., 2006. Development of Strategies for Improved Filtering and Fitting of SNMR Signals, MSc thesis, Technical University of Berlin. Walsh D. et al.  , 2013. A small-diameter NMR logging tool for groundwater investigations, Groundwater , 51( 6), 914– 926. https://doi.org/10.1111/gwat.12024 Google Scholar CrossRef Search ADS   Walsh D.O., 2008. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigations, J. Appl. Geophys. , 66( 3–4), 140– 150. https://doi.org/10.1016/j.jappgeo.2008.03.006 Google Scholar CrossRef Search ADS   Wang Z.X., 2009. The key technology of magnetic resonance sounding instrument for groundwater investigation, DSc thesis , Jilin University. Wu N., Li Y., Tian Y.N., Zhong T., 2016. Trace-transform-based time-frequency filtering for seismic signal enhancement in Northeast China, C. R. Geosci. , 348( 5), 360– 367. https://doi.org/10.1016/j.crte.2016.02.001 Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations