Open AccessArticle

Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation

Xiaokang Yao

^1,2

Jianmin Zhang

^2,3,

Zhenyang Yu

^1,2,

Fa Zhao

^1,2 and

Yong Sun

^1,2,*

College of Instrumentation and Electrical Engineering, Jilin University, Changchun 130061, China

Key Laboratory of Geo-exploration Instruments, Ministry of Education of China, Changchun 130061, China

College of Geo-Exploration Science and Technology, Jilin University, Changchun 130026, China

Author to whom correspondence should be addressed.

Remote Sens. 2019, 11(15), 1829; https://doi.org/10.3390/rs11151829

Submission received: 26 June 2019 / Revised: 22 July 2019 / Accepted: 27 July 2019 / Published: 5 August 2019

(This article belongs to the Special Issue Recent Advances in Subsurface Sensing Technologies)

Download

Browse Figures

Versions Notes

Abstract

The magnetic resonance sounding (MRS) method is a non-invasive, efficient and advanced geophysical method for groundwater detection. However, the MRS signal received by the coil sensor is extremely susceptible to electromagnetic noise interference. In MRS data processing, random noise suppression of noisy MRS data is an important research aspect. We propose an approach for intensive sampling sparse reconstruction (ISSR) and kernel regression estimation (KRE) to suppress random noise. The approach is based on variable frequency sampling, numerical integration and statistical signal processing combined with kernel regression estimation. In order to realize the approach, we proposed three specific sparse reconstructions, namely rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction. To solve the distortion of peaks and valleys after sparse reconstruction, we introduced the KRE to deal with the processed data by the ISSR. Further, the simulation and field experiments demonstrate that the ISSR-KRE approach is a feasible and effective way to suppress random noise. Besides, we find that rectangular sparse reconstruction and trapezoidal sparse reconstruction are superior to Simpson sparse reconstruction in terms of noise suppression effect, and sampling frequency is positively correlated with signal-to-noise improvement ratio (SNIR). In one case of field experiment, the standard deviation of noisy MRS data was reduced from 1200.80 nV to 570.01 nV by the ISSR-KRE approach. The proposed approach provides theoretical support for random noise suppression and contributes to the development of MRS instrument with low power consumption and high efficiency. In the future, we will integrate the approach into MRS instrument and attempt to utilize them to eliminate harmonic noise from power line.

Keywords:

magnetic resonance sounding (MRS); random noise suppression; sampling record method; signal estimation; intensive sampling; sparse reconstruction; kernel regression estimation

Graphical Abstract

1. Introduction

Magnetic resonance sounding (MRS) has attracted much attention in recent years due to accuracy properties of MRS in groundwater detection, which offers important economic benefits. MRS is a relatively efficient and non-invasive geophysical approach for groundwater surveying that utilizes the principle of magnetic resonance phenomenon in hydrogen atoms comprised of water in the geomagnetic field [1,2,3]. The development of nuclear magnetic resonance imaging in the biomedical field has promoted the MRS application in water resources exploration [4,5]. So far, the MRS method and instrument have been further researched and developed in many countries such as France, the United States, Germany, China and other countries [6,7,8]. The MRS method has been applied to hydrological investigation, observation of permafrost thaw, monitoring the cavity in glacier, water-induced disaster detection, and infiltrating water surveys [9,10,11,12,13]. However, when the MRS instrument is exploited for groundwater exploration, the received signal is particularly weak and easily degraded in ambient electromagnetic interference [14]. Furthermore, obtaining the MRS signal with nanovolt-level amplitude is a challenging area in regards to strong electromagnetic noise [15].

Due to noise interference, noise suppression is still a fundamental issue in MRS data processing and there has been extensive research regarding denoising methods in the past decades. The MRS signal is mainly disturbed by three kinds of noise, which are spike noise, harmonic noise and random noise. Spike noise mainly originates from instantaneous discharge, such as lightning, electric fence, and electronic switch. With regard to spike noise suppression, statistical stacking [16] and model-based subtraction [17] are used to deal with it. Harmonic noise mostly stems from power lines and it may seriously damage the MRS signal. In terms of harmonic noise suppression, there are several classical methods, such as bandpass filtering, notch filtering, adaptive notch filter, and model-based elimination method [18]. Nevertheless, notch filtering and adaptive notch filtering may cause the distortion of the phase and amplitude of the MRS signal. The model-based elimination needs to accurately estimate basic frequency of power harmonic noise. In addition, the method of dealing with co-frequency harmonics was studied by Wang et al. [19]. Apart from the interference from spike noise and harmonic noise, the influence of random noise is more common to MRS data. Random noise originates not only from the ambient environment, but also from the internal circuit of the MRS instrument. It is more complex to handle random noise than to cope with harmonic noise and spike noise. If random noise cannot be processed cleanly, it is difficult to obtain an accurate MRS signal. Hence, random noise suppression becomes more and more significant in MRS data processing. For the sake of analyzing the distribution of random noise, Dalgaard et al. developed a noise collector [20]. In order to suppress random noise, a stacking procedure and statistical stacking were used to improve the signal to noise ratio of the MRS signal [16]. Moreover, a global nonlinear least squares parameter fitting scheme was adopted to reduce signal bias amidst strong noise [21], but this scheme demanded synchronization detection which made the control of MRS instrument more complex. In particular, there is another method of eliminating random noise, namely adaptive noise elimination of multichannel MRS instrument [22]. However, this multi-channel noise elimination requires multiple reference coils, which increases the difficulty of field construction. Ghanati et al. proposed a randomized version of the singular value decomposition to accelerate singular spectrum analysis [23]. Besides, time-frequency peak filtering (TFPF) was applied to eliminate random noise of MRS signal [24]. TFPF has achieved good effect in processing random noise, but the signal processed by the method is an envelope of the MRS signal. Although these methods have solved some of the random noise interference in MRS groundwater detection, these methods do not always guarantee to eliminate random noise cleanly. Moreover, there are several denoising methods pertaining to hardware. For examples, the eight shape loop antenna was utilized to overcome high electromagnetic noise [25]. However, the denoising antenna reduces the depth of MRS instrument detection, which is only suitable for shallow groundwater survey. Behroozmand et al. used a central loop method to improve the signal to noise ratio (SNR) of magnetic resonance sounding data [26]. Furthermore, these methods do not involve the random noise elimination associated with sampling aspect.

If performing variable sampling or intensive sampling, then a more accurate signal can be obtained, but intensive sampling generates a lot of data. In data compression, signal and image feature extraction processing of other fields, they provide us some inspiration for processing MRS signals, such as synthetic aperture radar (SAR) images [27], penetrating radar noise filter [28], image processing [29] and so on. In particular, data compression and high-frequency sampling in SAR equipment provide us with exploring direction to process noisy MRS signal. Then, we ponder over the combination of these aspects with MRS data processing and noise suppression. To improve the work efficiency of groundwater detection and reduce the energy consumption of instrument, we propose a novel approach of intensive sampling sparse reconstruction (ISSR) and kernel regression estimation (KRE) to suppress random noise and obtain the MRS full wave signal. Then, to achieve this approach, the ISSR is specifically studied in depth and three sparse reconstructions are proposed and analyzed. We furtherly introduce kernel regression estimation into dealing with the MRS signal processed by the ISSR method. At last, the numerical simulations and field experiments are carried out to verify these methods. One of the outstanding advantages of the research is that suppressing random noise to improve SNR is studied based on a sampling perspective. On the one hand, because the ISSR method reduces the number of stacks by increasing the sampling rate, it provides a theoretical basis for the development of the MRS instrument with low power consumption and high work efficiency. On the other hand, three sparse reconstruction methods provide theoretical guidance and specific methods to suppress random noise by way of variable frequency sampling.

The rest of this paper is structured as follows: Section 2 presents characteristics of MRS signal and classical stacking method. Section 3 illustrates intensive sampling sparse reconstruction and kernel regression estimation approach for suppressing random noise. In Section 4, we conduct numerical simulations in-depth to verify and compare the theoretical denoising effect of the methods in Section 2 and Section 3. In Section 5, we implement field experiments to verify denoising performances of the methods, which are conducted in a magnetic shielding room and in Shaoguo Town, Changchun City, China respectively. In Section 6, we discuss the methods and point out limitations and potential application associated with the approach. In the final section, we provide the conclusions.

2. MRS Signal Analysis and the Classical Stacking Method

2.1. MRS Signal Analysis

Utilizing the MRS instrument for on-site groundwater detection, the signal collected by the receiver from the receiving coils includes both the MRS signal and electromagnetic noise. To further study noise elimination, it is necessary to understand the generation of MRS signal and the characteristics of the signal. The emission current causes the hydrogen protons in groundwater to jump toward higher energy levels. After the emission current is removed, hydrogen protons relax laterally and decay, and electromagnetic wave with Larmor frequency are radiated outward. Following this, the receiving coil receives the electromagnetic wave containing the MRS signal with Larmor frequency [1,30], as shown in Figure 1.

Neglecting the interference of ambient electromagnetic noise, the MRS full-wave signal

U (t)

induced in the receiving coil is expressed as Equation (1) [18].

U (t) = U_{0} \cdot e^{(- t / T_{2}^{*})} \cdot \cos (2 π f_{L a r m o r} t + φ_{0})

(1)

In Equation (1),

f_{L a r m o r}

denotes the Larmor frequency,

U_{0}

t

T_{2}^{*}

and

φ_{0}

are the initial amplitude, the recording time, the observed relaxation time and the initial phase respectively. Their units are Hz, nV, ms, ms and degree respectively.

In general, electromagnetic noise is divided into three categories: Spike noise, power frequency harmonic noise and random noise. The focus of this study is the suppression of random noise. Here we use

r (t)

to represent random noise, then the received noisy signal

S (t)

is expressed as follows:

S (t) = U (t) + r (t)

(2)

where the units of

S (t)

U (t)

, and

r (t)

are all nV.

2.2. Classic Stacking Method

In groundwater detection, the MRS signal belongs to the weak signal with nano-volt level. In the implementation of classical stacking method, hydrogen nuclei in groundwater are repeatedly excited by dozens of times, and dozens of data recordings are repeatedly performed at every pulse moment. Then, the collected data is stacked and averaged to achieve signal-to-noise improvement ratio (SNIR) increasing. To comprehend the SNIR concept and to analyze the performance of the method in depth, we present two forms of SNIR (voltage SNIR, power SNIR) of suppressing Gaussian noise according to the literature [31].

{\begin{cases} S N I R_{V} = \sqrt{N_{S}} \\ S N I R_{p} = N_{S} \end{cases}

(3)

where

N_{S}

represents the number of superposition times in the classic stacking method. Note that if the MRS signal is interfered by Gaussian noise, then the performing stacking of

N_{S}

times can make power SNR (

S N R_{P}

) improve by

N_{S}

times. That is, the power SNIR (

S N I R_{P}

) is proportional to the stacking number. However, if the noise is non-Gaussian, then the noise suppression effect is reduced by adopting the classic stacking method. In other words, if the noise is other types, then performing stacking of

N_{S}

times can make

S N R_{P}

improve by less than

N_{S}

times.

Although classical stacking can suppress the electromagnetic noise interference, there are several shortcomings associated with the method. The control of receiver and transmitter of MRS instrument needs high precision of time synchronization. Due to the multiple currents emitted by the MRS transmitter, the classical stacking method causes low working efficiency and high-power consumption of MRS instrument. The classical stacking method cannot satisfy the demands of fast and mobile detection of MRS technology. In other words, fast and mobile MRS instrument cannot stay too long above a location in groundwater detection.

3. Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation

3.1. Overall Approach of ISSR-KRE for Suppressing Random Noise

In order to overcome the deficiencies highlighted in the previous section of classical stacking method, we propose the intensive sampling sparse reconstruction and kernel regression estimation (ISSR-KRE) approach for suppressing random noise. The ISSR-KRE approach includes intensive sampling sparse reconstruction theory and kernel regression estimation theory. This approach only requires one or two excitations for each pulse moment in suppressing random noise. After each emission current is over, intensive sampling of the MRS signal, sparse reconstruction, and kernel regression estimation are performed to achieve random noise suppression.

The basic feature of the ISSR-KRE method is to adopt high sampling frequency

f_{s_{H}}

for data acquisition. On the premise that the sampling frequency

f_{s_{H}}

satisfies the sampling theorem (

f_{s_{H}} \geq 2 f_{c}

f_{c}

represents the highest frequency to measure and its unit is Hz.), if the sampling frequency

f_{s_{H}}

is higher, then more points of noisy signal can be obtained.

Although high sampling frequency can obtain plenty of signal information, it also has the disadvantages of containing a large amount of data as well as much more time taken in signal processing. Hence, the ISSR-KRE approach is proposed to achieve small amount of final data and random noise suppression, and the approach is implemented by the steps, as shown in Figure 2. The ISSR-KRE approach mainly consists of the following steps: (i) Setting the desired SNR, (ii) only collecting noise recording, (iii) at the first time, the signal is collected and estimated, (iv) determining to increase the sampling frequency based on the desired SNR, and (v) at the next time, the signal is collected and estimated.

3.2. Intensive Sampling Sparse Reconstruction for Suppressing Random Noise

3.2.1. Basic Frequency of Sparse Reconstruction

In intensive sampling sparse reconstruction, an important operation is to determine the basic frequency

f_{p r o p}

of sparse reconstruction. To obtain MRS signal waveform with higher fidelity after ISSR, we define the basic frequency

f_{p r o p}

, and its unit is Hz.

f_{p r o p}

can be understood as the sampling frequency of the reconstructed signal, which satisfies the sampling theorem and is less than

f_{s_{H}}

, as formula (4).

f_{s_{H}} > f_{p r o p} > 2 f_{c}

(4)

According to SNR requirements,

f_{s_{H}}

is calculated after determining

f_{p r o p}

based on formula (4). In order to prevent the frequency shift of MRS signal after reconstruction,

f_{p r o p}

needs further optimization. The optimization of

f_{p r o p}

is as shown in Equation (5).

f_{p r o p} = \arg {g (f_{p r o p}) \in Q}

(5)

where

\arg {•}

is an operator that represents solving the parameters that satisfy the relationship. That is, Equation (5) means finding the parameter

f_{p r o p}

that satisfies the relationship

g (f_{p r o p}) \in Q

In Equation (5),

{\begin{cases} g (f_{p r o p}) = \frac{f_{s_{H}}}{f_{p r o p}} \\ Q = {q | q = 2 n + 2, n \in N^{*}} \\ f_{p r o p} = k f_{L a r m o r} \\ k > 2 \end{cases}

(6)

where

g (f_{p r o p})

is a function of the argument

f_{p r o p}

Q

represents a set of integers,

q

is an element in the set

Q

N^{*}

is the positive integer set,

n

is an element in the positive integer set

N^{*}

, and

k

is a real number greater than 2.

After MRS sampling verification, we find that if the sampling rate

f_{s_{H}}

is relatively large and if the basic frequency

f_{p r o p}

satisfies

f_{p r o p} \geq 12 f_{L a r m o r}

, the reconstructed signal has higher fidelity. However,

f_{p r o p}

cannot be very large. Because

f_{p r o p}

is very large to cause a reduction in SNIR and data compression rate.

3.2.2. ISSR Implementation Process

After determining

f_{s_{H}}

and

f_{p r o p}

, the specific steps of ISSR is proposed in Figure 3. The ISSR method mainly includes: (i) Intensive sampling, (ii) sampled data segments, (iii) estimating the time and amplitude of each segment data to form a new reconstructed point, and (iv) sparse reconstructed signal is formed from the new reconstructed points calculated for each segment.

In Figure 3, the used sampling frequency is sampling frequency

f_{s_{H}}

with high frequency when performing intensive sampling. The length of each segment is

f_{s_{H}} / f_{p r o p}

in the data segmentation operation. Estimating the time and magnitude of each segmentation data is the core of the ISSR method to form new estimates. In estimating the time and magnitude process, random noise is suppressed and SNR is improved. In addition, another important advantage is that the data is compressed that helps to reduce storage costs and quickly calculate the inversion result of groundwater. In sparse reconstruction, the amount of data can be reduced by

f_{s_{H}} / f_{p r o p}

times. We note that the choice of reconstruction method is directly related to the signal fidelity and SNR improvement of the MRS signal. In the next few sections, we study the specific reconstruction methods to demonstrate their performance. Finally, the reconstructed signal is the signal with a high SNR and a small amount of data.

In order to comprehend intensive sampling of the ISSR method, we give Figure 4. In the figure,

j

denotes the number of time subinterval,

a

and

b

indicate the beginning and end of the time subinterval, respectively. Furthermore, we find that the adjacent sampling interval is very small when performing intensive sampling.

After defining the specific steps of intensive sampling sparse reconstruction, choosing what kind of sparse reconstruction to process each data segment becomes another key issue. To solve this issue, we attempt to employ integration theory. In the mean-value theorem for integrals, if a function

f (x)

is continuous on

[a, b]

then there exists

c

[a, b]

such that:

\frac{1}{b - a} \int_{a}^{b} f (x) d x = f (c)

. In other words, we can utilize integral median

f (c)

instead of

f (x)

[a, b]

. Inspired by this theorem, we explore to adopt a numerical integration method to sparsely reconstruct the signal. In numerical integration, there are three commonly used numerical integration methods: Numerical integration of the rectangular method, numerical integration of the trapezoidal method, and numerical integration of the Simpson method [32]. Moreover, the order of three integration methods, according to integration precision from superior to secondary, is Simpson numerical integration, trapezoidal numerical integration, and rectangular numerical integration.

So far, there are two problems: Whether it is feasible to perform intensive sampling sparse reconstruction of a MRS signal with random noise utilizing these three methods and which method is the best in three methods to suppress random noise. In order to solve the two problems, rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction are studied respectively. The reason why they are named is so that they are derived based on numerical integration methods. Below we obtain and analyze three methods to process noisy signal with random noise based on SNIR perspective.

3.2.3. Rectangular Sparse Reconstruction

After intensive sampling, the noisy signal is divided into segments with length n per segment where n equals

f_{s_{H}} / f_{p r o p}

. Then, according to intensive sampling sparse reconstruction process in Figure 3, we establish the theory of rectangular sparse reconstruction. The details can be found in Appendix A. Further, the specific rectangular sparse reconstruction is selected to process the noisy MRS signal. Importantly, a specific indicator, namely the signal to noise improvement ratio for voltage, is derived for evaluating noise elimination effect of rectangular sparse reconstruction.

S N I R_{V (Re ct)} = \frac{S N R_{V O}}{S N R_{V I}} = \sqrt{\frac{f_{s_{H}}}{f_{p r o p}}}

(7)

where SNR_VI is SNR for voltage prior to sparse reconstruction, SNR_VO is SNR for voltage after sparse reconstruction,

S N I R_{V (Re ct)}

denotes signal to noise improvement ratio for voltage of rectangular sparse reconstruction, and the unit of

S N I R_{V (Re ct)}

is one.

Then SNIR for power is expressed as Equation (8).

S N I R_{P (Re ct)} = \frac{f_{s_{H}}}{f_{p r o p}}

(8)

where,

S N I R_{P (Rect)}

is the signal to noise improvement ratio for power, and the unit of

S N I R_{P (Rect)}

is one.

It can be seen from Equation (8) that

S N I R_{P (Re ct)}

is proportional to

f_{s_{H}}

and inversely proportional to

f_{p r o p}

. Therefore, SNR can be improved by increasing the sampling frequency

f_{s_{H}}

or by selecting a reasonably small

f_{p r o p}

to suppress Gaussian random noise.

3.2.4. Trapezoidal Sparse Reconstruction

After intensive sampling, the noisy signal is divided into segments with length n per segment where n equals

f_{s_{H}} / f_{p r o p}

. Then, according to intensive sampling sparse reconstruction process in Figure 3, we establish the theory of trapezoidal sparse reconstruction. The details can be found in Appendix A. Further, the specific trapezoidal sparse reconstruction is selected to process the noisy MRS signal. Importantly, a specific indicator, namely signal to noise improvement ratio for voltage SNIR_V, is derived for evaluating noise elimination effect of trapezoidal sparse reconstruction.

S N I R_{V (Trap)} \approx \frac{2 f_{s_{H}}^{}}{\sqrt{4 f_{s_{H}} f_{p r o p} - 2 f_{s_{H}}^{2}}}

(9)

Then SNIR for power is expressed as Equation (10).

S N I R_{P (Trap)} \approx \frac{2 f_{s_{H}}^{2}}{2 f_{s_{H}} f_{p r o p} - f_{s_{H}}^{2}}

(10)

It can be seen from Equation (10) that

S N I R_{P (Trap)}

is positively correlated with

f_{s_{H}}

and negatively correlated with

f_{p r o p}

. Therefore, SNR can be improved by increasing the sampling frequency

f_{s_{H}}

or by selecting a reasonably small

f_{p r o p}

to suppress Gaussian random noise.

3.2.5. Simpson Sparse Reconstruction

After intensive sampling, the noisy signal is divided into segments with length n per segment where n equals

f_{s_{H}} / f_{p r o p}

. Then, according to intensive sampling sparse reconstruction process in Figure 3, we establish the theory of Simpson sparse reconstruction. Appendix A provides the details. Further, the specific Simpson sparse reconstruction is selected to process the noisy MRS signal. Importantly, a specific indicator, namely signal to noise improvement ratio for voltage, is derived for evaluating noise elimination effect of Simpson sparse reconstruction.

S N I R_{V (Simp)} \approx \frac{3 f_{s_{H}}^{}}{\sqrt{10 f_{s_{H}} f_{p r o p} - 2 f_{s_{H}}^{2}}}

(11)

Then SNIR for power is expressed as Equation (12).

S N I R_{P (Simp)} \approx \frac{9 f_{s_{H}}^{2}}{10 f_{s_{H}} f_{p r o p} - 2 f_{s_{H}}^{2}}

(12)

It can be seen from Equation (12) that

S N I R_{P (Simp)}

is positively correlated with

f_{s_{H}}

and negatively correlated with

f_{p r o p}

. Therefore, SNR can be improved by increasing the sampling frequency

f_{s_{H}}

or by selecting a reasonably small

f_{p r o p}

to suppress Gaussian random noise.

Ignoring the approximate calculation process influence, comparing the three sparse reconstruction methods according to Equations (8), (10), and (12), we can get: When

f_{s_{H}} / f_{p r o p}

is small, the order of noise suppression effect from superior to secondary is trapezoidal sparse reconstruction, rectangular sparse reconstruction, and Simpson sparse reconstruction. In addition, when

f_{s_{H}} / f_{p r o p}

is large, the noise elimination effect of rectangular sparse reconstruction is close to that of trapezoidal sparse reconstruction, and the noise elimination effect of Simpson sparse reconstruction is the worst. Considering the effect of approximate calculation on the denoising result, the comparison of three sparse reconstructions will be made in the simulation section.

3.3. Kernel Regression Estimation for Suppressing Random Noise

After intensive sampling and sparse reconstruction, there are distortion problems in the peak and valley positions of the signal, and SNR can still be further improved. In order to reduce the signal distortion, one-dimensional kernel regression is used to process the sparse reconstructed signal. The measurement model of one-dimensional signal can be expressed as Equation (13).

y_{i} = f (x_{i}) + ε_{i}, i = 1, 2, \dots, p

(13)

where

y_{i}

is the actual measured value of the signal at the i-th sampling point

x_{i}

f (•)

is the regression function,

x_{i}

denotes the time of the i-th sampling point,

ε_{i}

is the noise value, and

p

is the number of points participating in the regression calculation.

Although the specific form of

f (x_{i})

is uncertain, we can estimate it by a N-order Taylor series. If x is a point in the field of the sample point x_i, then the Taylor series expansion of the regression function

f (x_{i})

at the time

x

can be obtained based on Taylor series expansion theory [29].

\begin{array}{l} f (x_{i}) \approx f (x) + f^{'} (x) (x_{i} - x) + \frac{1}{2!} f^{″} (x) {(x_{i} - x)}^{2} + \dots + \frac{1}{N!} f^{(N)} (x) {(x_{i} - x)}^{N} \\ = β_{0} + β_{1} (x_{i} - x) + β_{2} {(x_{i} - x)}^{2} + \dots + β_{N} {(x_{i} - x)}^{N} \end{array}

(14)

According to the Taylor series above, if

β_{i}

is known, the regression function

f (x_{i})

can be obtained. In order to minimize the regression error, we utilize the local weighted least squares method to optimize.

\min_{{β_{n}}} \sum_{i = 1}^{P} [y_{i} - β_{0} - β_{1} (x_{i} - x) - β_{2} {(x_{i} - x)}^{2} - \dots - β_{N} {(x_{i} - x)}^{N}] \frac{1}{h} K (\frac{x_{i} - x}{h})

(15)

In Equation (15),

K (•)

is a kernel function that controls

x_{i} - x

weights.

| x_{i} - x |

is smaller, then the weight is larger.

h

is the smoothing parameter, which affects the penalty strength of

| x_{i} - x |

. When the smoothing factor h is larger, the regression curve of the signal becomes smoother; otherwise the regression curve is not smooth. The form of

K (•)

is not unique as long as the following conditions are met:

\int_{R^{1}} τ K (τ) d τ = 0, \int_{R^{1}} τ^{2} K (τ) d τ = c

(16)

where c is a constant. Common kernel functions include linear kernel function, Gaussian kernel function, and polynomial kernel function. Since the choice of kernel has little influence on the estimation accuracy, this paper chooses Gaussian kernel function to estimate the signal processed by sparse reconstruction. The Gaussian kernel function is as follows.

K_{h} (τ) = \frac{1}{h} \cdot \frac{1}{\sqrt{2 π}} \cdot e^{- \frac{1}{2} {(\frac{τ}{h})}^{2}}

(17)

where

τ

means

x_{i} - x

. For the selection of different orders N, local estimation of noisy signal can be achieved. In nonparametric statistical theory, N=0, 1, 2 is widely used.

(i) When

N = 0

, the famous Nadaraya-Watson estimate (NWE) is obtained, and its form is as shown in Equation (18). NWE is a local constant estimator.

\hat{f} (x_{i}) = β_{0} = \frac{\sum_{i = 1}^{P} K_{h} (x_{i} - x) y_{i}}{\sum_{i = 1}^{P} K_{h} (x_{i} - x)}, K_{h} (τ) = \frac{1}{h} K (\frac{τ}{h})

(18)

(ii) When

N = 1

, we can get a local linear kernel regression estimator.

\hat{f} (x_{i}) = X_{1} A^{- 1} B

(19)

where

X_{1} = [\begin{matrix} 1 & x_{i} \end{matrix} - x]

A = [\begin{matrix} \sum_{i = 1}^{P} K_{h} (x_{i} - x) & \sum_{i = 1}^{P} (x_{i} - x) K_{h} (x_{i} - x) \\ \sum_{i = 1}^{P} (x_{i} - x) K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{2} K_{h} (x_{i} - x) \end{matrix}]

B = [\begin{matrix} \sum_{i = 1}^{P} y_{i} K_{h} (x_{i} - x) \\ \sum_{i = 1}^{P} (x_{i} - x) y_{i} K_{h} (x_{i} - x) \end{matrix}]

(iii) When

N = 2

, we can get a local quadratic kernel regression estimator.

\hat{f} (x_{i}) = X_{2} A^{- 1} B

(20)

where

X_{2} = [{\begin{matrix} 1 & x_{i} - x & (x_{i} - x) \end{matrix}}^{2}]

A = [\begin{matrix} \sum_{i = 1}^{P} K_{h} (x_{i} - x) & \sum_{i = 1}^{P} (x_{i} - x) K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{2} K_{h} (x_{i} - x) \\ \sum_{i = 1}^{P} (x_{i} - x) K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{2} K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{3} K_{h} (x_{i} - x) \\ \sum_{i = 1}^{P} {(x_{i} - x)}^{2} K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{3} K_{h} (x_{i} - x) & \sum_{i = 1}^{P} {(x_{i} - x)}^{4} K_{h} (x_{i} - x) \end{matrix}]

B = [\begin{matrix} \sum_{i = 1}^{P} y_{i} K_{h} (x_{i} - x) \\ \begin{array}{l} \sum_{i = 1}^{P} (x_{i} - x) y_{i} K_{h} (x_{i} - x) \\ \sum_{i = 1}^{P} {(x_{i} - x)}^{2} y_{i} K_{h} (x_{i} - x) \end{array} \end{matrix}]

In general, N is positively correlated with the variance of estimated values and N is negatively correlated with the deviation of estimated values.

4. Numerical Simulations

The method proposed in this paper will be verified by simulation experiments and field experiments. In order to clearly compare the noise suppression in the experiments, the SNR for power, the mean squared error (MSE) and the SNIR for power (Equations (8), (10) and (12)) are utilized as indicators for evaluating the method effectiveness. SNR and MSE are defined as follows:

S N R = 10 \log_{10} (\frac{P_{U}}{P_{r}}) = 10 \log_{10} (\frac{\sum_{i = 1}^{N} U_{i}^{2}}{\sum_{i = 1}^{N} r_{i}^{2}})

(21)

where

P_{U}

and

P_{r}

represent the signal power and the noise power, respectively;

U_{i}

and

r_{i}

represent the signal voltage and the noise voltage, respectively. The unit of SNR for power is dB.

M S E = \frac{\sum_{i = 1}^{N} {[S_{i} - U_{i}]}^{2}}{N}

(22)

where

S_{i}

denotes the noisy signal voltage,

N

is the total number of discrete signal points. The unit of MSE is nV².

It can be seen from Equation (21) and Equation (22) that the larger the SNR is, the smaller the MSE is. Then,

S (t)

is closer to

U (t)

, and the noise reduction effect becomes better.

4.1. Simulations of Intensive Sampling Sparse Reconstruction for Random Noise Suppression

4.1.1. Sparse Reconstruction Simulation

In the simulation experiment, the acquisition time length is 250 ms, and the initial amplitude of the MRS signal

U_{0}

= 40.46 nV, the Larmor frequency

f_{L a r m o r}

= 2341.7 Hz,

T_{2}^{*}

= 500 ms. The synthesized noisy signal is composed of MRS signal and Gauss white noise. When synthesizing a noisy MRS signal, random noise is determined to be added into MRS signal based on the value of signal-to-noise ratio. In simulation, the initial SNR of synthetic noisy signal is −20 dB.

In Figure 5, the left-hand column is the time-series. The blue lines display the time-series of the sampling frequency

f_{p r o p}

. The red lines display the time-series of rectangular sparse reconstruction by different sampling frequencies. The right-hand column is the spectra. The blue lines display the spectra of the time-series of the sampling frequency

f_{p r o p}

, and the red lines display the spectra of the time-series of rectangular sparse reconstruction by different sampling frequencies that correspond to the left-hand column.

According to determining the basic frequency method in Section 3.2.1, firstly setting

f_{p r o p}

> 12

f_{L a r m o r}

, we calculate the solutions of

f_{p r o p}

by Equation (5) and Equation (6), then we choose the solution

f_{p r o p}

= 24

f_{L a r m o r}

. Following this, we adopt the sampling frequency

f_{s_{H}}

= 8

f_{p r o p}

, 16

f_{p r o p}

, 32

f_{p r o p}

, 64

f_{p r o p}

to intensively sample the noisy MRS signal and to restructure the noisy signal by way of rectangular sparse reconstruction. After rectangular sparse reconstruction, the SNR of the restructured signal is –10.9537 dB, –7.8869 dB, –4.9097 dB, –1.9680 dB and SNIR is 7.7996, 15.8032, 31.3668, 61.7511, respectively. The time domain waveform and the frequency domain amplitude of the reconstructed signals are shown in Figure 5.

The curves in Figure 6 are statistically obtained from 20 experimental simulations. In each simulation, the SNR of synthesized noisy signal is −20 dB. In Figure 6, the first row shows the curves of power signal to noise ratio by logarithmic representation and the curves of power signal-to-noise ratio by non-logarithmic representation, respectively. The unit of SNR by logarithmic representation is dB. The unit of SNR by non-logarithmic representation is 1. The second row shows the curves of SNIR for power and MSE, respectively. The unit of SNIR for power is 1. The unit of MSE is nV². The blue lines display simulation experimental value, and the red lines display the theoretical value.

From Figure 6, SNR and SNIR are proportional to the sampling frequency

f_{s_{H}}

when

f_{p r o p}

is constant. MSE is inversely proportional to

f_{s_{H}}

. The experiments show that the sparse reconstruction of rectangular method is feasible to suppress random noise, and Equations (7) and (8) are verified.

In addition, the simulation processes of trapezoidal sparse reconstruction and Simpson sparse reconstruction are similar to simulation process of rectangular sparse reconstruction. Here we only give simulation results of trapezoidal sparse reconstruction and Simpson sparse reconstruction.

In trapezoidal sparse reconstruction, when the basic frequency

f_{p r o p}

is constant, SNR and SNIR are positively related to the sampling frequency

f_{s_{H}}

and MSE is negatively related to the sampling frequency

f_{s_{H}}

. The experiments show that trapezoidal sparse reconstruction is feasible to suppress random noise, and Equations (9) and (10) are verified.

In Simpson sparse reconstruction, when the basic frequency

f_{p r o p}

is constant, SNR and SNIR are positively related to the sampling frequency

f_{s_{H}}

and MSE is negatively related to the sampling frequency

f_{s_{H}}

. The experiments show that the sparse reconstruction of Simpson method is feasible to suppress random noise, and Equations (11) and (12) are verified.

4.1.2. Comparison of Three Sparse Reconstructions

This sub-section compares rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction, and then we analyze the advantages and disadvantages of three sparse reconstructions. The comparison indicators include SNR, SNIR, MSE, and waveform.

In this simulation experiment, we see that the initial amplitude of noiseless MRS signal, the Larmor frequency, the observed relaxation time and the length of recording time are equal to 40.46 nV, 2341.7 Hz, 500 ms and 250 ms respectively. The SNR of synthesized noisy signal is −20 dB. The same collected data is processed by three different sparse reconstructions, and the indicators are compared.

Comparison of SNR, SNIR and MSE

In Figure 7, we compare the performance of three sparse reconstruction methods (namely (rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction). Figure 7a is comparison of signal-to-noise ratio of noisy MRS signal processed by the three methods. It can be seen from the figure that the signal-to-noise ratio of the noisy MRS signal processed by rectangular sparse reconstruction is better than that by other two kinds of sparse reconstruction. Figure 7b is a comparison of the signal-to-noise improvement ratio of a noisy MRS signal processed by three methods. We can find that the signal-to-noise improvement ratio of the noisy MRS signal processed by the rectangular sparse reconstruction is the best. Figure 7c is a comparison of the mean squared error of the noisy MRS signal processed by three methods. It can be seen that the MSE value of the noisy MRS signal processed by Simpson sparse reconstruction is larger than that by rectangular sparse reconstruction and trapezoidal sparse reconstruction.

In Table 1,

S N R_{Re c t}

S N R_{T r a p}

and

S N R_{S i m p}

represent the SNR of rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction, respectively. According to Figure 7 and Table 1, we can conclude that, in this case of strong noise, the order of three sparse reconstructions in the light of denoising effect from superior to secondary is rectangular sparse reconstruction, trapezoidal sparse reconstruction, and Simpson sparse reconstruction. With increasing sampling frequency, the denoising effect of trapezoidal sparse reconstruction is close to that of rectangular sparse reconstruction.

Waveform Comparison

From Figure 7 and Figure 8, we can get that: Three sparse reconstructions can suppress random noise and reconstruct the signal waveform. The denoising effect of rectangular sparse reconstruction is better than that of trapezoidal sparse reconstruction. The denoising effect of trapezoidal sparse reconstruction is superior to that of Simpson sparse reconstruction. With increasing sampling frequency

f_{s_{H}}

, the denoising effect of trapezoidal sparse reconstruction approaches the denoising effect of rectangular sparse reconstruction.

Different $f_{p r o p}$ Effect on Sparse Reconstruction

Through the comparisons in Figure 7 and Figure 8, the rectangular sparse reconstruction has good suppression noise effect among the three reconstruction methods. Choosing

f_{p r o p}

12 f_{L a m o r}

24 f_{L a m o r}

36 f_{L a m o r}

, we analyze the influence of different

f_{p r o p}

on the rectangular sparse reconstruction.

It can be seen from Figure 9 that, when Equation (5) is satisfied and the sampling frequency is a constant frequency, if the

f_{p r o p}

is smaller, then the SNR and SNIR are larger and the MSE is smaller after sparse reconstruction. Therefore, satisfying the fidelity condition of sparse reconstruction, a relatively small

f_{p r o p}

should be selected.

4.2. Comparison of Intensive Sampling Sparse Reconstruction and Classical Stacking Method

4.2.1. Simulations of Classical Stacking for Random Noise Suppression

In the simulation, the sampling frequency is 25 kHz, the acquisition time length is 250 ms, the initial amplitude of noiseless MRS signal is

U_{0}

= 40.46 nV, the Larmor frequency is

f_{L a r m o r}

= 2341.7 Hz, and the observed relaxation time is

T_{2}^{*}

= 500 ms. The synthesized noisy signal is composed of MRS signal and Gaussian white noise. When synthesizing the noisy MRS signal, random noise is determined to be added to MRS signal based on the value of signal-to-noise ratio. In simulation, the initial SNR of synthetic noisy signal is −20 dB.

In Figure 10, the left-hand column is the time-series. The blue lines display the unstacking time-series. The red lines display the stacking time-series. The right-hand column is the spectra. The blue lines display the spectra of the unstacking time-series, and the red lines display the spectra of the stacking time-series that correspond to the left-hand column.

After suppressing random noise by 8, 16, 32 and 64 times stacking data, the SNR for power is –10.9502 dB, –7.8772 dB, –4.9493 dB, –1.9423 dB and the SNIR is 7.8387, 15.9964, 32.3185, 63.2067, respectively. The time domain waveform and the frequency domain amplitude of the stacking signals are shown in Figure 10.

In Figure 11, the first row shows the curves of power signal to noise ratio by logarithmic representation and the curves of power signal-to-noise ratio by non-logarithmic representation, respectively. The unit of SNR by logarithmic representation is dB. The unit of SNR by non-logarithmic representation is 1. The second row shows the curves of SNIR for power and MSE, respectively. The unit of SNIR for power is 1. The unit of MSE is nV². The blue lines display simulation experimental value, and the red lines display theoretical value of the power signal-to-noise ratio. SNR is a non-logarithmic representation of the power signal to noise ratio.

From Figure 11, it can be concluded that SNR and SNIR are proportional to the number of stacks and MSE is inversely proportional to the number of stacks.

4.2.2. Comparison of Intensive Sampling Sparse Reconstruction and Classical Stacking Method

Through the above simulation experiments, the rectangular sparse reconstruction has the good effect of noise suppression. Therefore, we compare the intensive sampling and sparse reconstruction of rectangular method with the classical stacking method.

From Table 2 and Table 3, it can be seen that the SNR obtained after 64 times stacks is approximately equal to that of sparse reconstruction utilizing sampling frequency

f_{s_{H}} = 64 f_{p r o p}

According to Figure 5, Figure 6, Figure 10 and Figure 11, Table 2 and Table 3, the efficiency of ISSR method is higher than that of classical stacking method based on comparing the classical stacking method with the ISSR method. Furthermore, the sparse reconstruction method has a good denoising effect. Therefore, we can increase the sampling frequency instead of multiple long-term acquisitions. Furthermore, the ISSR method can be used to improve efficiency and reduce the power consumption of the MRS instrument.

4.3. Simulation of Kernel Regression Estimation

To obtain a better denoising effect and verify the validity of kernel regression estimation, the sparse reconstructed data is subjected to kernel regression estimation processing. It should be noted that the choice of kernel has little influence on the estimation accuracy. The Gaussian kernel function is a commonly used kernel function, so this part selects the Gaussian kernel function. After the kernel function is selected, the estimation effect is affected by three factors: The regression order, the window size of kernel function, and the smoothing factor.

4.3.1. Window Size and Smoothing Factor Effect on the Estimation Result

Window Size Effect

In analyzing the effect of window length on the estimation results of three different order regressions, we set the smoothing factor as a constant. The value of smoothing factor h is reasonably chosen based on a time interval of the reconstructed signal series. If the smoothing factor h is larger, the regression curve of the signal becomes smoother; otherwise the regression curve is not smooth. Here we choose the smoothing factor h to be equal to

5 / f_{p r o p}

where

1 / f_{p r o p}

is a time interval of the reconstructed signal utilizing ISSR method. In the analysis, the data to be processed by kernel regression is the data reconstructed by rectangular sparse reconstruction. The SNR of the data to be processed by kernel regression is 7.0737 dB, and

M S E

= 101.5047

{nV}^{2}

As shown in Figure 12a, the Nadaraya–Watson estimator and the linear kernel regression estimator have a good estimation effect when the window size is 7. By these two kernel regression estimators processing, the SNR of data is 13.5430 dB. Moreover, the quadratic kernel regression estimator has a good estimation effect when the window size equals 19. By quadratic kernel regression estimator processing, the SNR of data is 15.1425 dB.

Further analysis shows that quadratic kernel regression estimator is better than linear kernel regression estimator and Nadaraya–Watson estimator when the processed data has a relatively high SNR.

Smoothing Factor Effect

With constant window size to analyze the effect of the smoothing factor

h

on the estimation results, the window size of the Nadaraya–Watson estimator, linear kernel regression estimator and quadratic kernel regression estimator is 7, 7 and 19 respectively. In the analysis, the data to be processed by kernel regression is the data reconstructed by rectangular sparse reconstruction. The SNR of the data to be processed utilizing kernel regression is SNR = 7.0737 dB, and MSE = 101.5047 nV².

From Figure 12b, when smoothing factor

h

equals

2 / f_{p r o p}

, Nadaraya–Watson estimator (N = 0) and linear kernel estimator (N = 1) have a good estimation effect. By these two kernel regression estimators processing, the SNR is 13.8916 dB, 13.8907 dB, respectively. For the local quadratic estimator (N = 2), it has good estimation effect when

h

equals

5 / f_{p r o p}

. The SNR of data is 15.1425 dB by quadratic kernel regression estimator.

4.3.2. Waveform Comparison of Kernel Regression Estimation

The data to be processed by kernel regression is the data reconstructed by rectangular sparse reconstruction, and its

S N R

= 7.0737 dB,

M S E

= 101.5047

{nV}^{2}

. In the case of a high SNR, the local quadratic estimator can obtain a higher SNR, therefore selecting a local quadratic estimator, we set its window size to 19 and its smoothing factor

h = 5 / f_{p r o p}

for simulation experiment comparison. The waveforms and frequency domain amplitudes of kernel regression estimation are shown in Figure 13.

From Figure 12 and Figure 13, we can get that: Performing a local quadratic estimator of kernel function regression can further improve the signal SNR, such as adopting

f_{s} = 512 f_{p r o p}

and by intensive sampling sparse reconstruction and local quadratic estimator. SNR increases from the initial −20 dB to 15.1425 dB, and the SNIR reaches 3267.8. From the low frequency band (1–500 Hz) and the higher frequency band (4500–5000 Hz) in the spectrogram, and the spectrum amplitude near the Larmor frequency, we can find that the local estimator of kernel function regression has a low-pass filtering effect. Increasing SNR by kernel regression estimation also causes a small attenuation of the signal amplitude. However, this attenuation effect of kernel regression estimation can be reduced by selecting appropriate estimation parameters (namely, the regression order, the window size of kernel function, and the smoothing factor).

5. Field Experiments

In order to further verify the feasibility and effectiveness of the intensive sampling sparse reconstruction and kernel regression estimation, we conducted two experiments. The first experiment was carried out in laboratory and the experimental device consisted of signal generator and signal acquisition device. The second experiment was conducted in the field and the experimental setup was MRS instrument with high sampling rate.

5.1. Laboratory Experiment

In the experiment, the arbitrary waveform generator (AWG) and the signal collector were connected by a BNC (Bayonet Neill-Concelman) cable, and the signal collector acquired the noisy signal generated by the AWG signal generator. The signal collector recorded the noisy MRS signal for 250 ms at sampling frequency

f_{s}

= 500 kHz, and the number of sampling points is 125,000 points. The initial amplitude of MRS signal is 5 mV and the Larmor frequency

f_{L a r m o r}

is 2341.7 Hz. Calculating the solution

f_{p r o p} = 13.3450 f_{L a r m o r}

according to Equation (5) and

f_{L a r m o r}

, then the obtained noisy data containing random noise is processed by rectangular sparse reconstruction, as shown in Figure 14.

From Figure 14, it is concluded that the reasonable

f_{p r o p}

is calculated and selected, and the intensive sampling sparse reconstruction of the rectangular method can suppress noise and improve SNR. After the rectangular sparse reconstruction, SNR can be further improved by kernel regression estimation as described in Section 4.3, which is not described here.

5.2. Processing Experiment of Noisy MRS Data

The experimental data collection was carried out utilizing the MSR instrument in Shaoguo Town, Changchun City, China, as shown in Figure 15. The local Larmor frequency

f_{L a r m o r}

is 2335 Hz. The single-turn transmitting coil and the single-turn receiving coil are 100 m × 100 m squares. The sampling frequency of MRS receiver is 50 kHz and a single signal acquisition time is 256 ms. The maximum transmitting peak current is 208.7556 A and the emission duration is 40 ms. The single noisy signal collected is firstly eliminated by spike noise and removed power frequency harmonics of dozens of frequencies by an adaptive notch filter.

After removing spike noise and power-frequency noise of 50 Hz multiples from the collected noisy data, there is still residual noise in data. The residual noise is not completely random noise, as shown in Figure 16. After removing the power frequency harmonics, the signal data is processed by rectangular sparse reconstruction to verify the noise elimination effect of the ISSR method, as shown in Figure 16.

In Figure 16, since the sampling frequency is only 50 kHz, the smaller

f_{p r o p}

= 3.569

f_{L a r m o r}

= 8333.3 Hz is calculated and selected to perform rectangular sparse reconstruction according to Equation (5).

From comparison of low frequency sampling result and the ISSR method results in Figure 16, it can be concluded that: Although the denoising ability of rectangular sparse reconstruction is weakened in the case that the noise is non-Gaussian noise, rectangular sparse reconstruction still has a certain effect of suppressing noise. In other words, the ISSR method can suppress Gaussian random noise and other types of random noise.

Furtherly, the subsequent processing is performed to verify the effect of the ISSR method and extract the MRS signal parameters. Due to low SNR of the data processed by rectangular sparse reconstruction in the experiment, local kernel regression estimation is not utilized in this case that SNR is low. In non-Gaussian random noise where the denoising effect is not excellent, we adopt a bandpass filter (bandpass range [2325 Hz 2345 Hz]) to further process the data. Finally, we process the MRS signal segment from 30 ms to 150 ms by the least squares linear fitting method to extract the initial amplitude

U_{0}

= 425.42 nV and the observed relaxation time

T_{2}^{*}

= 136.54 ms.

6. Discussion

We provide an approach to suppress random noise by utilizing high sampling frequency in the MRS instrument. In the approach, the basic frequency

f_{p r o p}

of sparse reconstruction is one of the ISSR keys. This basic frequency is accurately calculated according to Equation (5) to avoid frequency offset and distortion of MRS signal. In the case where the sampling frequency is constant, the magnitude of basic frequency

f_{p r o p}

is inversely related to the SNIR of signal in intensive sampling sparse reconstruction, as shown in Figure 9.

Three sparse reconstruction methods provide theoretical basis and specific methods for high frequency sampling to suppress random noise. These methods are feasible and have good effects according to Figure 5, Figure 6, Figure 14 and Figure 16. The rectangular sparse reconstruction and trapezoidal sparse reconstruction are better than Simpson sparse reconstruction in terms of the suppression effect of random noise (especially, in suppressing Gaussian random noise), from Figure 7a–c and Table 1 in Section 4 and Equations (7)–(12) in Section 3.2.3 for details. This order is the opposite of the order of the three numerical integration effects. This phenomenon of order inconsistency is caused by noise interference. By comparing the classical stacking method and the intensive sampling sparse reconstruction in the paper, we can increase sampling frequency instead of multiple long-term recordings of classical stacking method, and improve the efficiency of the instrument, seeing Figure 5, Figure 6, Figure 10 and Figure 11, Table 2 and Table 3. Besides, these sparse reconstruction methods lay a theoretical foundation for the compression of MRS data with high sampling frequency. The sparse reconstruction method has a good effect in suppressing random noise, but it also has limitations such as poor performance in suppressing harmonic noise and requiring high sampling frequency. In terms of sparse reconstruction limitation in Section 5.2, for example, when there is a large amount of non-Gaussian noise such as harmonic noise, the effect of sparse reconstruction method is weakened and it needs combination with other filtering methods.

The kernel regression estimation introduced in this paper provides a feasible method for suppressing noise in MRS data with high SNR, seeing Figure 13. To obtain better regression estimation results, we should pay attention to the selection of reasonable window size, order, and the smoothing factor according to Figure 12. In addition, the regression estimation does not work well in the case of large noise.

The intensive sampling sparse reconstruction approach has fewer recording times, less transmitting energy consumption and higher efficiency in terms of instrument transmission and reception. Furthermore, the approach provides theoretical guidance and specific methods for high or variable frequency sampling to suppress random noise. At the same time, the approach gives a research direction for the development of MRS instruments with low power consumption and high work efficiency. In addition, the intensive sampling sparse reconstruction approach can also be applied to suppressing random noise in other geophysical signal processing or other signal processing areas. In the future, we will endeavour to overcome the limitations of ISSR-KRE and to handle other types of noise (i.e., harmonic noise, spike noise) in combination with other methods.

7. Conclusions

In order to better suppress random noise in MRS data, we propose an approach of intensive sampling sparse reconstruction and kernel regression estimation to improve SNR in this paper. The approach is based on variable frequency sampling, numerical integration and statistical signal processing combined with kernel regression estimation. Furthermore, we propose three specific sparse reconstructions (i.e., rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction) to realize intensive sampling sparse reconstruction. By means of simulative experiments, we find that all three sparse reconstructions can suppress random noise, and rectangular sparse reconstruction and trapezoidal sparse reconstruction are better than Simpson sparse reconstruction in terms of the suppression effect of random noise. At the same time, we verify the positive correlation of SNIR and sampling frequency. Besides, compared with the classical stacking method, the sparse reconstruction method has high efficiency and good denoising effect. In a simulation example, the ISSR is utilized to process noisy MRS data, and SNR is improved by 27.0737 dB. Further, in simulation of kernel regression estimation, we find that the suppression effect of kernel regression is related to Taylor order, window size and smoothing factor. Hence, we should choose kernel regression estimation with appropriate parameters to reduce random noise. In coping with the data processed by the ISSR, SNR is increased by 8.0688 dB adopting the KRE. In other words, the ISSR-KRE approach improved SNR by 35.1425 dB. Furthermore, in both the laboratory experiment and the field experiment of the MRS instrument, we verify the feasibility and effectiveness of the ISSR-KRE approach to suppress random noise by processing the data sampled at different sampling frequencies, such as 500 kHz and 50 kHz. For example, in processing noisy MRS data, the standard deviation of the data is reduced from 1200.80 nV to 570.01 nV utilizing the ISSR. The proposed approach and specific methods provide theoretical support for random noise suppression. In addition, the approach contributes to sampling method and the development of MRS instrument with low power consumption and high efficiency. In the future, we will integrate the approach and the methods into the hardware system of MRS instrument and attempt to utilize these methods to eliminate power frequency noise.

Author Contributions

Conceptualization, X.Y.; methodology, X.Y. and J.Z.; data curation Z.Y. and Y.S.; software, Y.S. and F.Z.; formal analysis, Z.Y., J.Z, and Y.S.; writing—original draft, X.Y. and Z.Y.; writing—review and editing, J.Z., X.Y. and Y.S.; project administration, Y.S.; funding acquisition, X.Y.

Funding

This work was supported by the National Natural Science Foundation for Excellent Young Scholars of China under Grant No. 41722405, the National Key Foundation for Exploring Scientific Instrument of China under Grant No. 2017YFC0804105, and the Key Research and Development Program of Jilin Province of China under Grant No. 20180201017GX and under Grant No. 20160414002GH.

Acknowledgments

Thanks are due to Tingting Lin, Ling Wan, Yang Zhang and Yujing Yang for valuable discussions and suggestions, and to Deli Sun, Chunjuan He and Jinbao Zhu for assistance with the experiments. We would like to thank the editor and reviewers for their reviews that improved the content of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Sparse Reconstruction

(1) Rectangular sparse reconstruction

After intensive sampling, the whole data is divided into segments with length n per segment. In other words, each data segment has

n

small segments, as in Equation (A1).

n = \frac{f_{s_{H}}}{f_{p r o p}}

(A1)

As shown in Figure 4, in the interval

[t_{a j}, t_{b j}]

the data is equally divided. In other words, the interval

[t_{a j}, t_{b j}]

is divided into

n

subintervals with equal length by the points

t_{a j} = t_{0 j}, t_{1 j}, t_{2 j}, \dots, t_{n j} = t_{b j}

, and the length of each subinterval is

Δ t = (t_{b j} - t_{a j}) / n

. Here,

j

denotes the number of interval.

i

represents the number of one point in interval. a and b describe the starting point and ending point of an interval, respectively.

In a subinterval

[t_{(i - 1) j}, t_{i j}]

of interval

[t_{a j}, t_{b j}]

, the area of approximate rectangle surrounded by the data curve and subinterval

[t_{(i - 1) j}, t_{i j}]

[(t_{b j} - t_{a j}) / n] \cdot U_{i j}

. According to the rectangular integral method and Equation (A1), Equation (A2) is obtained.

\int_{t_{a j}}^{t_{b j}} U (t) d t \approx \frac{t_{b j} - t_{a j}}{n} (U_{1 j} + U_{2 j} + \dots + U_{n j})

(A2)

Further, when

(t_{b j} - t_{a j}) < 1 / f_{p r o p}

, in order to estimate the signal amplitude

{\hat{U}}_{j}

[t_{a j}, t_{b j}]

, we consider the integrated area of the signal

U (t)

in time

[t_{a j}, t_{b j}]

as a rectangle area with length

{\hat{U}}_{j}

and width

t_{b j} - t_{a j}

, so

{\hat{U}}_{j}

can be estimated.

{\hat{U}}_{j} = \frac{\int_{t_{a j}}^{t_{b j}} U (t) d t}{t_{b j} - t_{a j}} \approx \frac{\sum_{i = 1}^{n} U_{i j}}{n}

(A3)

This rectangular estimation is now utilized to study random suppression and the improvement of SNR. Since white noise is the common and important random noise, we assume that the standard deviation of the zero-mean Gaussian distribution random noise

r_{i j}

σ_{r}

, i.e.,

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

, then for a single point sampling

S_{i j} = U_{i j} + r_{i j}

where the MRS signal is

U_{i j}

, SNR for voltage at

t_{i j}

is Equation (A4).

S N R_{V i j} = U_{i j} / σ_{r}

(A4)

Then, in the time domain

[t_{a j}, t_{b j}]

the average SNR for voltage of the noisy signal is Equation (A5). It is SNR prior to sparse reconstruction.

S N R_{V I j} = \frac{S N R_{V 1 j} + S N R_{V 2 j} + \dots + S N R_{V n j}}{n} = \frac{\sum_{i = 1}^{n} U_{i j}}{n σ_{r}}

(A5)

The signal estimated by the rectangular method is

{\hat{S}}_{j}

[t_{a j}, t_{b j}]

, as Equation (A6).

{\hat{S}}_{j} \approx \frac{S_{1 j} + S_{2 j} + \dots + S_{n j}}{n} = \frac{\sum_{i = 1}^{n} S_{i j}}{n}

(A6)

Since

S_{i j} = U_{i j} + r_{i j}

, the above equation can be expressed as Equation (A7).

{\hat{S}}_{j} \approx \frac{\sum_{i = 1}^{n} S_{i j}}{n} = \frac{\sum_{i = 1}^{n} U_{i j} + \sum_{i = 1}^{n} r_{i j}}{n} = {\hat{U}}_{j} + \frac{1}{n} \sum_{i = 1}^{n} r_{i j}

(A7)

According to

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

and the statistical laws that linear combinations of a finite number of mutually independent Gaussian random variables still obey Gaussian distributions [33], after using sparse reconstruction

S N R_{V O j}

can be derived and obtained, as Equation (A8).

S N R_{V O j} = \sqrt{n} \cdot \frac{\hat{U_{j}}}{σ_{r}}

(A8)

[t_{a j}, t_{b j}]

SNIR for voltage can be derived and obtained according to Equations (A1), (A5) and (A8).

S N I R_{V} = S N R_{V O j} / S N R_{V I j} = \sqrt{\frac{f_{s_{H}}}{f_{p r o p}}}

(A9)

Then SNIR for power is expressed as

S N I R_{P}

, Equation (A10).

S N I R_{P (Re ct)} = \frac{f_{s_{H}}}{f_{p r o p}}

(A10)

(2) Trapezoidal sparse reconstruction

First, the trapezoidal sparse reconstruction is used to study the signal without noise. As shown in Figure 4, using the points

t_{a j} = t_{0 j}, t_{1 j}, t_{2 j}, \dots, t_{n j} = t_{b j}

to divide

[t_{a j}, t_{b j}]

into

n

equal parts, the equidistant points

t_{k j} = t_{a} + k Δ t, (k = 0, 1, \dots, n)

are selected and each interval length is

Δ t = (t_{b j} - t_{a j}) / n

n = \frac{f_{s_{H}}}{f_{p r o p}}

(A11)

In each subinterval

[t_{k j}, t_{(k + 1) j}]

, the trapezoidal formula is used, and the compound trapezoidal formula (A12) is obtained.

\int_{t_{a j}}^{t_{b j}} U (t) d t = \sum_{k = 0}^{n - 1} \int_{t_{k j}}^{t_{(k + 1) j}} U (t) d t \approx \frac{Δ t}{2} \sum_{k = 0}^{n - 1} [U (t_{k}) + U (t_{k + 1})]

(A12)

Considering the integrated area of the signal

U (t)

in time

[t_{a j}, t_{b j}]

as a rectangle area with length

{\hat{U}}_{j}

and width

t_{b j} - t_{a j}

, so

{\hat{U}}_{j}

can be estimated.

{\hat{U}}_{j} = \frac{\int_{t_{a j}}^{t_{b j}} U (t) d t}{t_{b j} - t_{a j}} \approx \frac{\sum_{k = 0}^{n - 1} [U_{k j} + U_{(k + 1) j}]}{2 n}

(A13)

This trapezoidal method estimation is now utilized to study the improvement of SNR in the case where the acquired signal contains Gaussian random noise. Assuming that the standard deviation of the zero-mean Gaussian distribution random noise

r_{i j}

σ_{r}

, i.e.,

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

, then for a single point sampling

S_{i j} = U_{i j} + r_{i j}

where the MRS signal is

U_{i j}

, SNR for voltage at

t_{i j}

is Equation (A14).

S N R_{V i j} = U_{i j} / σ_{r}

(A14)

Then, the average SNR for voltage of the signal is Equation (A15) in the time domain

[t_{a j}, t_{b j}]

. It is SNR prior to sparse reconstruction.

S N R_{V I j} = \frac{S N R_{V 1 j} + S N R_{V 2 j} + \dots + S N R_{V n j}}{n} = \frac{\sum_{i = 1}^{n} U_{i j}}{n σ_{r}}

(A15)

The signal estimated by the trapezoidal method is

{\hat{S}}_{j}

[t_{a j}, t_{b j}]

, as Equation (A16).

{\hat{S}}_{j} \approx \frac{\sum_{k = 0}^{n - 1} [S_{k j} + S_{(k + 1) j}]}{2 n} = \frac{\sum_{k = 0}^{n - 1} [U_{k j} + U_{(k + 1) j}]}{2 n} + \frac{\sum_{k = 0}^{n - 1} [r_{k j} + r_{(k + 1) j}]}{2 n}

(A16)

According to

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

and the statistical laws that linear combinations of a finite number of mutually independent Gaussian random variables still obey Gaussian distributions [33], Equation (A17) is derived.

\frac{1}{2 n} \sum_{k = 0}^{n - 1} [r_{k j} + r_{(k + 1) j}] \sim N (0, \frac{2 n - 1}{2 n^{2}} σ_{r}^{2})

(A17)

After sparse reconstruction,

S N R_{V O j}

can be derived and obtained, as in Equation (A18).

S N R_{V O j} = \frac{\frac{\sum_{k = 0}^{n - 1} [U_{k j} + U_{(k + 1) j}]}{2 n}}{\sqrt{\frac{2 n - 1}{2 n^{2}} σ_{r}^{2}}} = \frac{\sum_{k = 0}^{n - 1} [U_{k j} + U_{(k + 1) j}]}{σ_{r} \sqrt{4 n - 2}} = \frac{2 n \hat{U_{j}} + U_{0 j} - U_{n j}}{σ_{r} \sqrt{4 n - 2}}

(A18)

[t_{a j}, t_{b j}]

SNIR for voltage can be derived and obtained according to Equations (A15) and (A18).

S N I R_{V} = \frac{S N R_{V O j}}{S N R_{V I j}} = \frac{(2 n \hat{U_{j}} + U_{0 j} - U_{n j})}{\hat{U_{j}} \sqrt{4 n - 2}} = \frac{2 n}{\sqrt{4 n - 2}} + \frac{U_{0 j} - U_{n j}}{\hat{U_{j}} \sqrt{4 n - 2}}

(A19)

When n is relatively large and

| U_{0 j} - U_{n j} | ≪ | \hat{U_{j}} |

, then the term

\frac{U_{0 j} - U_{n j}}{\hat{U_{j}} \sqrt{4 n - 2}}

can be tended to zero. In detail, if the sampling is very dense in the interval

[t_{a j}, t_{b j}]

, then the value of n is very large. If

| U_{0 j} - U_{n j} |

is far less than

| \hat{U_{j}} |

, then

\frac{U_{0 j} - U_{n j}}{{\hat{U}}_{j}}

can be tended to zero. Further, when n is relatively large and

| U_{0 j} - U_{n j} | ≪ | \hat{U_{j}} |

, we can get

\frac{2 n}{\sqrt{4 n - 2}} / \frac{| U_{0 j} - U_{n j} |}{| \hat{U_{j}} | \sqrt{4 n - 2}} ≫ 1

, and the formula (A20) is obtained by Equation (A11).

S N I R_{V} \approx \frac{2 n}{\sqrt{4 n - 2}} = \frac{2 f_{s_{H}}}{\sqrt{4 f_{s_{h}} f_{p r o p} - 2 f_{p r o p}^{2}}}

(A20)

Then SNIR for power is expressed as

S N I R_{P}

, Equation (A21).

S N I R_{P (Trap)} \approx \frac{2 f_{s_{H}}^{2}}{2 f_{s_{H}} f_{p r o p} - f_{s_{H}}^{2}}

(A21)

(3) Simpson sparse reconstruction

Firstly, the trapezoidal sparse reconstruction is utilized to study the signal without noise. As shown in Figure 4, using the points

t_{a j} = t_{0 j}, t_{1 j}, t_{2 j}, \dots, t_{n j} = t_{b j}

to divide

[t_{a j}, t_{b j}]

into

n

equal parts, each interval length is

Δ t = \frac{t_{b j} - t_{a j}}{n}

n = \frac{f_{s_{H}}}{f_{p r o p}}

(A22)

According to the Newton–Cotes formula (an interpolation type of integral formula), the Simpson formula in the interval

[t_{(i - 1) j}, t_{(i + 1) j}]

is derived.

\int_{t_{(i - 1) j}}^{t_{(i + 1) j}} U (t) d t \approx \frac{Δ t}{3} [U (t_{(i - 1) j}) + 4 U (t_{i j}) + U (t_{(i + 1) j})]

(A23)

Further, the compound Simpson formula in the interval

[t_{a j}, t_{b j}]

is obtained.

\int_{t_{b j}}^{t_{b j}} U (t) d t \approx \frac{Δ t}{3} [U_{0 j} + U_{n j} + 4 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} U_{i j} + 2 \sum_{\begin{array}{l} i = 2 k \\ k = 1, 2 \dots \end{array}}^{n - 2} U_{i j}]

(A24)

Considering the integrated area of the signal

U (t)

in time

[t_{a j}, t_{b j}]

as a rectangle area with length

{\hat{U}}_{j}

and width

t_{b j} - t_{a j}

, so

{\hat{U}}_{j}

can be estimated.

{\hat{U}}_{j} = \frac{\int_{t_{a j}}^{t_{b j}} U (t) d t}{t_{b j} - t_{a j}} \approx \frac{1}{3 n} [U_{0 j} + U_{n j} + 4 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} U_{i j} + 2 \sum_{\begin{array}{l} i = 2 k \\ k = 1, 2 \dots \end{array}}^{n - 2} U_{i j}]

(A25)

Further, we utilize the Simpson method estimation to study the improvement of SNR in the case where the acquired noisy signal contains Gaussian random noise. Assuming that the standard deviation of the zero-mean Gaussian distribution white noise

r_{i j}

σ_{r}

, i.e.,

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

, then for a single point sampling

S_{i j} = U_{i j} + r_{i j}

where the MRS signal is

U_{i j}

, SNR for voltage at

t_{i j}

is Equation (A26).

S N R_{V i j} = U_{i j} / σ_{r}

(A26)

Then, the average SNR for voltage of the noisy signal is Equation (A27) in the time domain

[t_{a j}, t_{b j}]

. It is SNR prior to sparse reconstruction.

S N R_{V I j} = \frac{S N R_{V 1 j} + S N R_{V 2 j} + \dots + S N R_{V n j}}{n} = \frac{\sum_{i = 1}^{n} U_{i j}}{n σ_{r}}

(A27)

The signal estimated by the Simpson method is

{\hat{S}}_{j}

[t_{a j}, t_{b j}]

, as Equation (A28).

{\hat{S}}_{j} \approx \frac{1}{3 n} [U_{0 j} + U_{n j} + 4 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} U_{i j} + 2 \sum_{\begin{array}{l} i = 2 k \\ k = 1, 2 \dots \end{array}}^{n - 2} U_{i j}] + \frac{1}{3 n} [r_{0 j} + r_{n j} + 4 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} r_{i j} + 2 \sum_{\begin{array}{l} i = 2 k \\ k = 1, 2 \dots \end{array}}^{n - 2} r_{i j}]

(A28)

According to

r_{i j} \sim N (0, σ_{r_{i j}}^{2})

and the statistical laws that linear combinations of a finite number of mutually independent Gaussian random variables still obey Gaussian distributions [33], Equation (A29) is derived.

\frac{1}{3 n} [r_{0 j} + r_{n j} + 4 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} r_{i j} + 2 \sum_{\begin{array}{l} i = 2 k \\ k = 1, 2 \dots \end{array}}^{n - 2} r_{i j}] \sim N (0, \frac{10 n - 2}{9 n^{2}} σ_{r}^{2})

(A29)

After Simpson method sparse reconstruction

S N R_{V O j}

can be derived and obtained, as Equation (A30).

S N R_{V O j} = \frac{2 n \hat{U_{j}} + 2 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} U_{i j} + U_{0 j} - U_{n j}}{σ_{r} \sqrt{10 n - 2}}

(A30)

[t_{a j}, t_{b j}]

SNIR for voltage can be derived and obtained according to Equations (A27) and (A30).

S N I R_{V} = \frac{2 n}{\sqrt{10 n - 2}} + \frac{2 \sum_{\begin{array}{l} i = 1 + 2 k \\ k = 0, 1, \dots \end{array}}^{n - 1} U_{i j}}{\hat{U_{j}} \sqrt{10 n - 2}} + \frac{U_{0 j} - U_{n j}}{\hat{U_{j}} \sqrt{10 n - 2}}

(A31)

When

n

is relatively large and

| U_{0 j} - U_{n j} | ≪ | \hat{U_{j}} |

, formula (A32) is obtained by Equations (A22) and (A31).

S N I R_{V} \approx \frac{3 f_{s_{H}}}{\sqrt{10 f_{s_{H}} f_{p r o p} - 2 f_{_{p r o p}}^{2}}}

(A32)

Then SNIP for power is expressed as

S N I R_{P}

, Equation (A33).

S N I R_{P (Simp)} \approx \frac{9 f_{s_{H}}^{2}}{10 f_{s_{H}} f_{p r o p} - 2 f_{s_{H}}^{2}}

(A33)

References

Hertrich, M.; Braun, M.; Gunther, T.; Green, A.G.; Yaramanci, U. Surface nuclear magnetic resonance tomography. IEEE Trans. Geosci. Remote Sens. 2007, 45, 3752–3759. [Google Scholar] [CrossRef]
Legchenko, A.; Valla, P. Removal of power-line harmonics from proton magnetic resonance measurements. J. Appl. Geophys. 2003, 53, 10–120. [Google Scholar] [CrossRef]
Liu, L.; Grombacher, D.; Auken, E.; Larsen, J.J. Removal of Co-Frequency Powerline Harmonics from Multichannel Surface NMR Data. IEEE Geosci. Remote Sens. Lett. 2018, 15, 53–57. [Google Scholar] [CrossRef]
Chen, H.M.; Wang, H.C.; Chai, J.W.; Chen, C.C.C.; Xue, B.; Wang, L.; Yu, C.; Wang, Y.; Song, M.; Chang, C.I. A Hyperspectral Imaging Approach to White Matter Hyperintensities Detection in Brain Magnetic Resonance Images. Remote Sens. 2017, 9, 1174. [Google Scholar] [CrossRef]
Powers, J.M.; Ioachim, G.; Stroman, P.W. Ten Key Insights into the Use of Spinal Cord fMRI. Brain Sci. 2018, 8, 173. [Google Scholar] [CrossRef] [PubMed]
Legchenko, A.; Descloitres, M.; Bost, A.; Ruiz, L.; Reddy, M.; Girard, J.F.; Sekhar, M.; Kumar, M.S.M.; Braun, J.J. Resolution of MRS applied to the characterization of hard-rock aquifers. Groundwater 2006, 44, 547–554. [Google Scholar] [CrossRef] [PubMed]
Walsh, D.O. Multi-channel surface NMR instrumentation and software for 1D/2D groundwater investigations. J. Appl. Geophys. 2008, 66, 140–150. [Google Scholar] [CrossRef]
Qin, S.; Ma, Z.; Jiang, C.; Lin, J.; Xue, Y.; Shang, X.; Li, Z. Response Characteristics and Experimental Study of Underground Magnetic Resonance Sounding Using a Small-Coil Sensor. Sensors 2017, 17, 2127. [Google Scholar] [CrossRef]
Valois, R.; Vouillamoz, J.M.; Lun, S.; Arnout, L. Mapping groundwater reserves in northwestern Cambodia with the combined use of data from lithologs and time-domain-electromagnetic and magnetic-resonance soundings. Hydrogeol. J. 2018, 26, 1187–1200. [Google Scholar] [CrossRef]
Parsekian, A.D.; Creighton, A.L.; Jones, B.M.; Arp, C.D. Surface nuclear magnetic resonance observations of permafrost thaw below floating, bedfast and transitional ice lakes. Geophysics 2019, 84, EN33–EN45. [Google Scholar] [CrossRef]
Garambois, S.; Legchenko, A.; Vincent, C.; Thibert, E. Ground-penetrating radar and surface nuclear magnetic resonance monitoring of an englacial water-filled cavity in the polythermal glacier of Tete Rousse. Geophysics 2016, 81, WA131–WA146. [Google Scholar] [CrossRef]
Shang, X.; Jiang, C.; Ma, Z.; Qin, S. Combined System of Magnetic Resonance Sounding and Time-Domain Electromagnetic Method for Water-Induced Disaster Detection in Tunnels. Sensors 2018, 18, 3508. [Google Scholar] [CrossRef] [PubMed]
Falzone, S.; Keating, K. Algorithms for removing surface water signals from surface nuclear magnetic resonance infiltration surveys. Geophysics 2016, 81, WB97–WB107. [Google Scholar] [CrossRef]
Ghanati, R.; Hafizi, M.K.; Fallahsafari, M. Surface nuclear magnetic resonance signals recovery by integration of a non-linear decomposition method with statistical analysis. Geophys. Prospect. 2016, 64, 489–504. [Google Scholar] [CrossRef]
Liu, L.; Grombacher, D.; Auken, E.; Larsen, J.J. Complex envelope retrieval for surface nuclear magnetic resonance data using spectral analysis. Geophys. J. Int. 2019, 217, 894–905. [Google Scholar] [CrossRef]
Jiang, C.; Lin, J.; Duan, Q.; Sun, S.; Tian, B. Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements. Near Surf. Geophys. 2011, 9, 459–468. [Google Scholar] [CrossRef]
Larsen, J.J. Model-based subtraction of spikes from surface nuclear magnetic resonance data. Geophysics 2016, 81, WB1–WB8. [Google Scholar] [CrossRef] [Green Version]
Larsen, J.; Dalgaard, E.; Auken, E. Noise cancelling of MRS signals combining model-based removal of powerline harmonics and multichannel Wiener filtering. Geophys. J. Int. 2014, 196, 828–836. [Google Scholar] [CrossRef]
Wang, Q.; Jiang, C.; Mueller-Petke, M. An alternative approach to handling co-frequency harmonics in surface nuclear magnetic resonance data. Geophys. J. Int. 2018, 215, 1962–1973. [Google Scholar] [CrossRef]
Dalgaard, E.; Christiansen, P.; Larsen, J.J.; Auken, E. A temporal and spatial analysis of anthropogenic noise sources affecting SNMR. J. Appl. Geophys. 2014, 110, 34–42. [Google Scholar] [CrossRef]
Legchenko, A.; Valla, P. Processing of surface proton magnetic resonance signals using non-linear fitting. J. Appl. Geophys. 1998, 39, 77–83. [Google Scholar] [CrossRef]
Dalgaard, E.; Auken, E.; Larsen, J.J. Adaptive noise cancelling of multichannel magnetic resonance sounding signals. Geophys. J. Int. 2012, 191, 88–100. [Google Scholar] [CrossRef]
Ghanati, R.; Hafizi, M.K.; Mahmoudvand, R.; Fallahsafari, M. Filtering and parameter estimation of surface-NMR data using singular spectrum analysis. J. Appl. Geophys. 2016, 130, 118–130. [Google Scholar] [CrossRef]
Lin, T.; Zhang, Y.; Yi, X.; Fan, T.; Wan, L. Time-frequency peak filtering for random noise attenuation of magnetic resonance sounding signal. Geophys. J. Int. 2018, 213, 727–738. [Google Scholar] [CrossRef]
Trushkin, D.V.; Shushakov, O.A.; Legchenko, A. The potential of a noise-reducing antenna for surface NMR groundwater surveys in the Earth’s magnetic field. Geophys. Prospect. 1994, 42, 855–862. [Google Scholar] [CrossRef]
Behroozmand, A.A.; Auken, E.; Fiandaca, G.; Rejkjaer, S. Increasing the resolution and the signal-to-noise ratio of magnetic resonance sounding data using a central loop configuration. Geophys. J. Int. 2016, 205, 243–256. [Google Scholar] [CrossRef] [Green Version]
Karine, A.; Toumi, A.; Khenchaf, A.; El Hassouni, M. Radar Target Recognition Using Salient Keypoint Descriptors and Multitask Sparse Representation. Remote Sens. 2018, 10, 843. [Google Scholar] [CrossRef]
Zhang, J.; Zeng, Z.; Zhang, L.; Lu, Q.; Wang, K. Application of Mathematical Morphological Filtering to Improve the Resolution of Chang’E-3 Lunar Penetrating Radar Data. Remote Sens. 2019, 11, 524. [Google Scholar] [CrossRef]
Takeda, H.; Farsiu, S.; Milanfar, P. Kernel Regression for Image Processing and Reconstruction. IEEE Trans. Image Process. 2007, 16, 349–366. [Google Scholar] [CrossRef] [Green Version]
Legchenko, A. Magnetic Resonance Imaging for Groundwater; ISTE Ltd: London, UK, 2013. [Google Scholar]
Legchenko, A.; Valla, P. A review of the basic principles for proton magnetic resonance sounding measurements. J. Appl. Geophys. 2002, 50, 3–19. [Google Scholar] [CrossRef]
Li, Q.; Wang, N.; Yi, D. Numerical Analysis, 5th ed.; Tsinghua University Press: Beijing, China, 2008; pp. 97–135. [Google Scholar]
Sheng, Z.; Xie, S.; Pan, C. Probability Theory and Mathematical Statistics, 4th ed.; Higher Education Press: Beijing, China, 2008; pp. 76–83. [Google Scholar]

Figure 1. Schematic diagram for detection of magnetic resonance sounding (MRS) signal.

Figure 2. Flowchart for implementation of the intensive sampling sparse reconstruction and kernel regression estimation (ISSR-KRE) approach.

Figure 3. Schematic diagram of intensive sampling sparse reconstruction.

Figure 4. A data fragment of intensive sampling in

[t_{a j}, t_{b j}]

Figure 4. A data fragment of intensive sampling in

[t_{a j}, t_{b j}]

Figure 5. Simulation case: Comparison of suppressing random noise effect of rectangular sparse reconstruction by different sampling frequencies

f_{s_{H}}

. Time-series and spectra of the sampling frequencies 8

f_{p r o p}

, 16

f_{p r o p}

, 32

f_{p r o p}

, and 64

f_{p r o p}

, respectively.

Figure 5. Simulation case: Comparison of suppressing random noise effect of rectangular sparse reconstruction by different sampling frequencies

f_{s_{H}}

. Time-series and spectra of the sampling frequencies 8

f_{p r o p}

, 16

f_{p r o p}

, 32

f_{p r o p}

, and 64

f_{p r o p}

, respectively.

Figure 6. The curves of signal to noise ratio (SNR), signal-to-noise improvement ratio (SNIR), mean squared error (MSE) with increasing sampling frequency by rectangular sparse reconstruction.

Figure 7. Comparison of three sparse reconstruction methods (rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction): (a) SNR comparison, (b) SNIR comparison, (c) MSE comparison.

Figure 8. Comparison of rectangular sparse reconstruction, trapezoidal sparse reconstruction and Simpson sparse reconstruction for processing signals. (a) Time domain, (b) frequency domain.

Figure 9. Different

f_{p r o p}

effect on sparse reconstruction. The first row shows the comparison of SNR in different

f_{p r o p}

cases, respectively. The second row shows the comparison of SNIR and MSE in different

f_{p r o p}

cases, respectively.

Figure 9. Different

f_{p r o p}

effect on sparse reconstruction. The first row shows the comparison of SNR in different

f_{p r o p}

cases, respectively. The second row shows the comparison of SNIR and MSE in different

f_{p r o p}

cases, respectively.

Figure 10. Simulation case: Time-series and spectra of 8, 16, 32, 64 stacks at one pulse moment, respectively.

Figure 11. The curves of SNR, SNIR and MSE with increasing stacking times.

Figure 12. The effect of kernel regression parameters on estimation results. (a) The effect of window size. (b) The effect of the smoothing factor h.

Figure 13. Waveform comparison of the data processed by the local quadratic estimator (window size 19 and smoothing factor

h = 5 / f_{p r o p}

). (a) Comparison in time domain. (b) Comparison in frequency domain. The green lines display the data by the rectangular sparse reconstruction method. The red lines display the ideal signal. The blue lines display the data by the rectangular sparse reconstruction method and the local quadratic kernel regression estimation. In other words, the blue lines display the ISSR-KRE data.

Figure 13. Waveform comparison of the data processed by the local quadratic estimator (window size 19 and smoothing factor

h = 5 / f_{p r o p}

Figure 14. Comparison of traditional low frequency sampling and the ISSR method. The gray lines display traditional low frequency sampling recording. The blue lines display the data processed by the ISSR method. (a) Comparison in time domain. (b) Comparison in frequency domain.

Figure 15. Field experiment location and the MRS instrument.

Figure 16. Comparison of low frequency sampling result and the ISSR method results. The green lines display low frequency sampling recording (its sampling frequency is 8330 Hz) and the standard deviation is 1200.80 nV. The blue lines display the data by the ISSR method (its sampling frequency is 50 kHz,

f_{p r o p}

= 8333.3 Hz) and the standard deviation is 570.01 nV. (a) Comparison in time domain. (b) Comparison in frequency domain.

f_{p r o p}

= 8333.3 Hz) and the standard deviation is 570.01 nV. (a) Comparison in time domain. (b) Comparison in frequency domain.

Table 1. SNR comparison of three ISSR methods in the experiment. The estimated parameters,

S N R_{Re c t}

S N R_{T r a p}

and

S N R_{S i m p}

are reported based on the synthetic MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a m o r}

= 2341.7 Hz,

f_{p r o p}

= 24

f_{L a m o r}

and initial SNR = −20 dB.

Table 1. SNR comparison of three ISSR methods in the experiment. The estimated parameters,

S N R_{Re c t}

S N R_{T r a p}

and

S N R_{S i m p}

are reported based on the synthetic MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a m o r}

= 2341.7 Hz,

f_{p r o p}

= 24

f_{L a m o r}

and initial SNR = −20 dB.

$f_{s_{H}} (Hz)$	$8 f_{p r o p}$	$16 f_{p r o p}$	$32 f_{p r o p}$	$64 f_{p r o p}$	$128 f_{p r o p}$	$256 f_{p r o p}$	$512 f_{p r o p}$
$S N R_{Re c t} - S N R_{T r a p} (d B)$	0.2519	0.0999	0.0610	0.0487	0.0253	0.0032	0.0027
$S N R_{Re c t} - S N R_{S i m p} (d B)$	0.3419	0.3650	0.3974	0.4625	0.4689	0.4814	0.4970

Table 2. The classical stacking method. The estimated parameters, SNR and SNIR are reported based on the synthetic noisy MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a r m o r}

= 2341.7 Hz, and initial SNR = −20 dB.

Table 2. The classical stacking method. The estimated parameters, SNR and SNIR are reported based on the synthetic noisy MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a r m o r}

= 2341.7 Hz, and initial SNR = −20 dB.

$N_{S}$	1	8	16	32	64
$S N R (d B)$	−20	−10.9502	−7.8772	−4.9493	−1.9423
$S N I R$	1	7.8387	15.9964	32.3185	63.2067

Table 3. The intensive sampling rectangular sparse reconstruction. The estimated parameters, SNR and SNIR are reported based on the synthetic MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a r m o r}

= 2341.7 Hz,

f_{p r o p}

= 24

f_{L a r m o r}

and SNR = −20 dB.

Table 3. The intensive sampling rectangular sparse reconstruction. The estimated parameters, SNR and SNIR are reported based on the synthetic MRS signal with

U_{0}^{}

= 40.46 nV,

T_{2}^{*}

= 500 ms,

f_{L a r m o r}

= 2341.7 Hz,

f_{p r o p}

= 24

f_{L a r m o r}

and SNR = −20 dB.

$f_{s_{H}} (Hz)$	$f_{p r o p}$	$8 f_{p r o p}$	$16 f_{p r o p}$	$32 f_{p r o p}$	$64 f_{p r o p}$
$S N R (dB)$	−20	−10.9537	−7.8869	−4.9097	−1.9680
$S N I R$	1	7.7996	15.8032	31.3668	61.7511

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Yao, X.; Zhang, J.; Yu, Z.; Zhao, F.; Sun, Y. Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation. Remote Sens. 2019, 11, 1829. https://doi.org/10.3390/rs11151829

AMA Style

Yao X, Zhang J, Yu Z, Zhao F, Sun Y. Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation. Remote Sensing. 2019; 11(15):1829. https://doi.org/10.3390/rs11151829

Chicago/Turabian Style

Yao, Xiaokang, Jianmin Zhang, Zhenyang Yu, Fa Zhao, and Yong Sun. 2019. "Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation" Remote Sensing 11, no. 15: 1829. https://doi.org/10.3390/rs11151829

APA Style

Yao, X., Zhang, J., Yu, Z., Zhao, F., & Sun, Y. (2019). Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation. Remote Sensing, 11(15), 1829. https://doi.org/10.3390/rs11151829

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

Article Menu

Random Noise Suppression of Magnetic Resonance Sounding Data with Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation

Abstract

1. Introduction

2. MRS Signal Analysis and the Classical Stacking Method

2.1. MRS Signal Analysis

2.2. Classic Stacking Method

3. Intensive Sampling Sparse Reconstruction and Kernel Regression Estimation

3.1. Overall Approach of ISSR-KRE for Suppressing Random Noise

3.2. Intensive Sampling Sparse Reconstruction for Suppressing Random Noise

3.2.1. Basic Frequency of Sparse Reconstruction

3.2.2. ISSR Implementation Process

3.2.3. Rectangular Sparse Reconstruction

3.2.4. Trapezoidal Sparse Reconstruction

3.2.5. Simpson Sparse Reconstruction

3.3. Kernel Regression Estimation for Suppressing Random Noise

4. Numerical Simulations

4.1. Simulations of Intensive Sampling Sparse Reconstruction for Random Noise Suppression

4.1.1. Sparse Reconstruction Simulation

4.1.2. Comparison of Three Sparse Reconstructions

Comparison of SNR, SNIR and MSE

Waveform Comparison

Different f p r o p Effect on Sparse Reconstruction

4.2. Comparison of Intensive Sampling Sparse Reconstruction and Classical Stacking Method

4.2.1. Simulations of Classical Stacking for Random Noise Suppression

4.2.2. Comparison of Intensive Sampling Sparse Reconstruction and Classical Stacking Method

4.3. Simulation of Kernel Regression Estimation

4.3.1. Window Size and Smoothing Factor Effect on the Estimation Result

Window Size Effect

Smoothing Factor Effect

4.3.2. Waveform Comparison of Kernel Regression Estimation

5. Field Experiments

5.1. Laboratory Experiment

5.2. Processing Experiment of Noisy MRS Data

6. Discussion

7. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

Appendix A. Sparse Reconstruction

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Different $f_{p r o p}$ Effect on Sparse Reconstruction