Research on NMR Noise Reduction Method Based on Improved CEEMD

Low-field NMR technology has been widely used in many fields, among which the most representative method is the use of CPMG sequence for T2 measurement. Important parameters such as permeability, saturation and fluid type of porous media can be obtained by T2 spectrum inversion calculation of NMR echo sequence signal. However, T2 spectrum inversion calculation is easily affected by noise signals. Improving the signal-to-noise ratio (SNR) of echo signal is one of the key problems with the application of low-field NMR technology. Therefore, in order to obtain more accurate analysis results, it is necessary to develop corresponding noise reduction methods. The past researches on NMR signal noise reduction mainly focus on echo string. In this paper, a complete NMR signal model is constructed, and the problem of the echo and echo string de-noising is discussed. FIR filter and orthogonal modulation filter (OMF) are studied for noise reduction of echo. To solve the problem of echo string de-noising, an improved complementary ensemble empirical mode decomposition (CEEMD) threshold de-noising method is proposed and compared with the heuristic threshold de-noising method based on sym4 wavelet and the previous empirical mode decomposition (EMD) threshold de-noising method. Based on the construction model and practical verification, it is proved that the proposed method can invert the more accurate T2 spectrum while improving the signal-to-noise ratio. The research results of this paper provide a strong support for the application of low-field NMR technology in the environment of strong noise.


I. INTRODUCTION
Low-field NMR technology characterized by relaxation measurement has been widely used in petroleum exploration, biomedicine, quality control and other fields, among which the most representative method is the use of CPMG sequence [1] for T2 measurement. In order to obtain more accurate T2 spectrum, it is necessary to improve the signalto-noise ratio of NMR logging signals.
The existing NMR noise reduction methods generally adopt the different distribution of signal and noise in time domain, frequency domain or time-frequency domain to reduce noise. Due to their advantages in NMR signal processing, these noise reduction methods are taken seriously [2]- [5]. In particular, the wavelet noise reduction method, which can maintain the peak signal while noise reduction features, has become a current research hotspot.
The associate editor coordinating the review of this manuscript and approving it for publication was Fang Yang . In recent years, Huang proposed a new time-frequency analysis method -empirical mode decomposition (EMD). These methods have been widely used in signal noise reduction in various fields [6]- [9]. EMD is an adaptive decomposition method. But it has the problem of modal aliasing. With the study of EMD, Huang proposed the integrated empirical mode decomposition method (EEMD) to improve the EMD method by adding white noise to better suppress the phenomenon of modal aliasing. However, white noise added in EEMD cannot be completely neutralized, so the researchers propose complementary ensemble empirical mode decomposition (CEEMD). CEEMD does not only guarantee the decomposition effect, but also reduces the reconstruction error caused by white noise.
The existing NMR signal de-noising methods generally assume that the noise is white noise subject to Gaussian distribution, but in fact it is much more complicated. For the de-noising of NMR echo string, many scholars adopt wavelet threshold de-noising method, and put forward many improved methods based on it [14], [15]. However, the wavelet threshold de-noising method is based on the underlying white noise model, which cannot eliminate the colored noise. Furthermore, the appropriate wavelet base and decomposition layers number should be determined before wavelet decomposition. This greatly increases the complexity and difficulty of the application and affects the application effect. CEEMD is a decomposition based on raw data. Therefore, this paper proposes a CEEMD method based on improved threshold value for echo string noise reduction. Fractional Gaussian noise (FGn) was used to model the noise of the echo string. The Hurst parameter [16] of FGn was estimated to build the threshold function to effectively remove all kinds of noises. In this paper, the proposed method and the existing methods are simulated and compared to verify the advantages of this method. The measured data are processed to verify the effectiveness of the method.
In this paper, in order to achieve the optimal de-noising effect, noise reduction methods are studied for both the spin echo signal and the echo string signal. The rest of the paper is organized as follows: the second section determines the overall noise reduction scheme and introduces the noise reduction methods involved in it in detail. In the third section, numerical simulation and experiment are carried out to evaluate the noise reduction effect on the change of SNR and the comparison of inversion effect. The fourth section summarizes the work of this paper.

II. METHOD
The SNR of the echo series signal directly affects the accuracy of relaxation time of the inversion. The overall scheme flow chart for noise reduction is determined as shown in Figure 1. For noise reduction of spin echos, we adopt FIR filtering and Orthogonal modulation filtering-OMF (abbreviated as FIR-OMF), respectively extracting peak information about spin echos after each filtering to form echo string. In this paper, CEEMD improved threshold method is adopted to further filter and reduce noise of echo string signal. Compared with wavelet threshold noise reduction and EMD threshold noise reduction, the advantages of the proposed method are verified.

A. FIR-OMF SPIN ECHO NOISE REDUCTION
The noise reduction flow chart of each spin echo signal is shown in Figure 2, where I is the same phase component, Q is the orthogonal component, and 2(Q 2 n + I 2 n ) 1/2 is the envelope of each echo. After noise reduction by FIR-OMF, we can get the echo string by taking the peak value of the envelope. The FIR-OMF spin echo noise reduction process is described as follows: (1) Each spin echo signal is filtered by FIR bandpass filter. The noise component is filtered out. Then the signal multiplies by sine component and cosine component of the same frequency, respectively, to get the orthogonal component q n and the same direction component i n .
(2) Pass q n and i n through the FIR low-pass filter, filter the frequency multiplication component, get Q n and I n .
(3) According to formula 2(Q 2 n + I 2 n ) 1/2 , the envelope of each filtered spin echo signal is obtained.

B. NOISE REDUCTION OF ECHO STRING
For noise reduction of echo string, CEEMD threshold method is adopted in this paper. The core of CEEMD threshold noise reduction method is to construct an appropriate threshold according to the noise energy estimated by each IMF.

1) BRIEF REVIEW OF CEEMD
Essentially, CEEMD method is the average result of multiple EMD decomposition of adaptive additive Gaussian white noise. By adding small white noise to make the signal balanced, the phenomenon of mode aliasing is effectively solved. Meanwhile, the real signal is maintained based on the zero mean property of Gaussian white noise, which is a big improvement on the traditional EMD analysis method.
The procedure of CEEMD can be described as follows.
(1) First, white noise is added to the original signal to obtain the first EMD component of the data with noise. Each time a different noise is added to the EMD decomposition. The decomposition is repeated many times. The average value of the first IMF component obtained from the decomposition of EMD is used as the first IMF component of the decomposition of CEEMD, i.e., Here, VOLUME 8, 2020 generating the i th IMF component, and N is the number of decomposition.
(2) Second, calculate the first residual signal r 1 (t), i.e., ..N ) as the new signal for EMD decomposition. Similar to the step (1), the average value of the first IMF component obtained from the decomposition of EMD in this step is used as the second IMF component of the CEEMD decomposition, i.e., (4) Finally, repeat the above steps; thus, we can obtain the (j+1) th IMF component C j+1 (t), i.e., R(t) is the last residual function. Through the above process, CEEMD decomposition and reconstruction are realized.

2) FRACTIONAL GAUSSIAN NOISE MODEL
As previously mentioned, the core of CEEMD threshold noise reduction method is to construct an appropriate threshold according to the noise energy estimated by each IMF. Therefore, Fractional Gaussian noise (FGn) is introduced in this paper.
FGn is a generalization of discrete Gaussian white noise [17]- [21]. It is a universal model of uniform broadband noise without the main frequency band and is essentially a discrete time process [22]. The statistical characteristic of FGn is only determined by its second order structure, which only depends on the real value parameter H . FGn is generally defined as a zero-mean Gaussian stationary process. Its autocorrelation sequence is shown in Eq (5): where, τ is the lag length and σ 2 is the variance. The proposed CEEMD threshold de-noising method uses FGn to model the noise of echo string signal. As the core of FGn, the parameter H needs to be estimated to build the threshold function.
In this paper, periodic graph method was used to estimate H . Periodic graph is a basic numerical method of power spectral density (PSD) calculation of time series. For a discrete time series, its periodic diagram is shown in Eq (6): The PSD of FGn is obtained by applying the discrete Fourier transform into the auto-correlation function (5) as follows: where f is frequency and C is constant. If H is not equal to 0.5, the PSD of FGn is approximate to: Convert Eq (8) into logarithmic form: As shown in Eq (9), the slope of the logarithmic periodic graph corresponding to the logarithmic frequency is ρ = 1 − 2H . We can get the logarithmic periodic graph of the echo string signal to be processed, and then use the least square fitting method to fit it to get the fitting line. According to the slope of the fitting line, we get the H parameter we need.

3) SELECTION OF THRESHOLD FUNCTION
In this study, the EMD method was employed to decompose Fractal Gaussian noise and analyze the energy distribution of the noise signal on each IMF. The noise energy distribution of follow-up IMF can be estimated by using the noise energy of the first and the second component of the IMF(imf (1,2) (t)) [23]. Therefore, the threshold noise reduction method based on CEEMD firstly needs to calculate the noise energy of imf (1,2) (t). The estimation formula is as follows: The improved threshold is selected for filtering, and the specific expression of the threshold is as follows: Th j = σ j (2 ln N ) 1/2 / ln(2j + 1) j = 1, 2 (11) where N is the length of the processed signal. The relation between power spectral density (PSD) of noise distribution of each IMF and coefficient j is [16]: where, ρ H can be solved by the following quadratic equation: where H is the value of Hurst index between 0 and 1, and when the fractal Gaussian noise is Gaussian white noise H = 0.5 [24]. Then we can integrate the power spectral density of each IMF over the frequency spectrum to obtain the variation on noise power with the coefficient on each IMF: Since the value range of H is in (0,1), the existing relation (j > j ≥ 2) proves that the noise energy at each IMF decreases with the increase of the coefficient. At the same time, it can be found that the noise energy at any IMF can be calculated by using the energy at imf 2 (t).
Therefore, the threshold value on any IMF of j > j ≥ 2 can be calculated by the following formula: The improved CEEMD de-noising threshold is In the process of threshold noise reduction, CEEMD threshold noise reduction method is different from wavelet threshold noise reduction method. In general, hard threshold or soft threshold equation is not adopted.
For the application of threshold noise reduction to NMR signal, the proposed CEEMD threshold de-noising method is more inclined to hard threshold. However, direct threshold causes discontinuity in the thresholded IMFs. We adopt the more sophisticated interval threshold instead. According to the nature of IMF, all IMF quantities are composed of oscillations passing zero. A point on the j th IMF (imf (j) (t)) is denoted by d j (k) and all zero points of imf (j) (t) is denoted by Z m j . The threshold filtering method can be expressed as; when Z m j is in the scope of [Z m j , Z m+1 j ), if there is any |d j (k)| greater than Th j , we keep all the Z m j in [Z m j , Z m+1 j ). The interval threshold equation is shown as: where d j (Z m j ) indicates the samples of the interval Z m j , and r m j represents the single extremum of the interval [25]. The improved CEEMD threshold de-noising method is shown as Figure 3. Based on the analysis of the proposed method, its time complexity is O(n 3 ).

III. RESULTS
Numerical simulation and experiments are carried out in this section. In order to evaluate the effect of noise reduction, it is analyzed from the following two aspects: first, the change of signal-to-noise ratio; the second is the comparison of inversion effects.

A. NUMERICAL SIMULATIONS
According to the overall noise reduction scheme, as shown in Figure 4, we first built the simulation model.

1) ESTABLISHMENT OF SIMULATION MODEL
In the process of NMR logging, the stratigraphic information is very complex. For example, the pore size presents a distribution trend, and the liquid composition is also very complex. In the process of NMR logging, an echo string can be obtained by employing CPMG pulse echo attenuation. The string is not affected by a single T2 value, but with a certain functional relationship with different T2 values. The functional relationship is expressed as After the establishment of T2 spectrum model with twopeak, an attenuation echo string curve C1 was obtained. The process of constructing simulation signal is as follows: (1) A single echo signal is simulated based on the function (2) Extend a single echo signal to form 300 echoes.
(3) Multiply the result from step 1 by C1 to get the desired signal.
(4) Based on step 2, the simulation signal is obtained by adding noise.
Based on above, the proposed de-noising method includes the following steps: (1) Establish the simulation signal model.
(2) FIR and OMF are used to reduce the noise of each echo respectively. The peak value of each echo is extracted to form an echo string.
(3) Implement CEEMD for the echo string signal to obtain IMFs. VOLUME 8, 2020 (4) Estimate the standard deviation of the first two IMF components and the variance of the other IMF components by Eqs. (10) and (16). Calculate threshold values by Eqs. (11) and (17).
In FIR filtering, high order and high resolution are the cause of well effect for noise reduction. However, excessively high order may lead to the reduction of signal amplitude, which resulted in the lost of useful information. Therefore, in this model, the bandpass filter order is 32, and the lowpass filter order is 64.
This method is used for noise reduction of NMR signal. The effect is shown in Figure 5- Figure 13. Figure 5 (a) is the model diagram of the pure NMR signal. Figure 5 (b) is the amplified spin echo. Figure 6 (a) shows the noise-contaminated NMR signal with a signal-to-noise ratio of −33.3191dB. Figure 6 (b) shows the amplified spin echo. Figure 7 is the echo string signal formed by extracting the peak value of each spin echo signal after FIR-OMF filtering. Figure 8 is the fitting signal of logarithmic periodic graph of    the threshold noise reduction algorithm based on EMD and the improved threshold noise reduction method based on CEEMD are compared. The results are shown in Figure 9. In order to show that the method proposed in this paper can adapt to the environment of low SNR and is superior to other methods, according to the above analysis, we add noise with an SNR of −51.7415 to the pure signal for simulation. Similarly, Figure 10 and Figure 11 are the noise reduction effect diagrams with an initial SNR of −51.7415dB. Figure 13 is the reconstructed signal. Figure 12 is the fitting signal of the logarithmic periodic diagram of the signal shown in Figure 11, and H is 0.5474.

2) SNR ANALYSIS
In order to make a more accurate interpretation of the effect after noise reduction, the signal-to-noise ratio (SNR) is now  taken as the evaluation standard. In this paper, the following formula is used to calculate signal SNR: where, s(n) refers to the ideal noise-free signal, s(n) is the reconstructed signal, and N is the length of the signal.

VOLUME 8, 2020
The SNR changes after noise reduction of −33.3191dB signal and −51.7415dB signal were calculated respectively. The results are shown in Table 1. It can be seen that the FIR-OMF method improves SNR and achieves certain de-noising effect. For the de-noising of echo string, the three methods mentioned in this paper can obviously improve the signal-tonoise ratio. Under the background of strong noise, CEEMD method is obviously better than wavelet in improving SNR. This shows that the proposed method is suitable for strong noise reduction.

3) INVERSION ANALYSIS
The purpose of signal noise reduction is to obtain relatively accurate T2 spectrum. So it is very important to carry out inversion analysis on the reconstructed signal. In this paper, SVD inversion method is used to carry out inversion analysis of the various de-noising methods mentioned in this paper, and is compared with the ideal T2 spectrum. Figure 14 is an ideal T2 spectrum. In order to facilitate analysis, the T2 spectrum obtained by inversion of various methods is finally displayed on the graph, as shown in Figure 15 and Figure 16. From the inversion results, it can be seen that the FIR-OMF algorithm can improve the accuracy of T2 spectrum with the noise reduction of the spin echo. For noise reduction of echo string, the improved threshold noise reduction algorithm  based on CEEMD proposed in this paper is superior to other algorithms.

B. EXPERIMENTAL VALIDATION
The effectiveness of our proposed method is verified by practical experiments.
The experimental samples are shown in Figure 17, which are tight sandstone samples. Two-dimensional T1-T2 spectrum experiments were carried out on these samples with the 2MHz benchtop NMR spectrometer from Magritek, to obtain two-dimensional data, which was equivalent to measuring multiple groups of T2 in the process of longitudinal magnetization vector attenuation. The first three groups of them are taken for experimental analysis. Since the experimental instrument has embedded noise reduction algorithm, the original three sets of data are taken as ideal data. In addition, the three sets of data were processed by the method mentioned in this paper.   Signal to noise ratio and T2 spectrum inversion results are also used to verify the results. Table 2 shows the SNR comparison processed by the three methods. It can be seen   that the algorithm proposed in this paper improves the SNR compared with the other two methods. Figure 18- Figure 20 are the T2 spectrum inversion results after processing the three groups of data respectively. It can be seen that the inversion T2 spectrum obtained by the proposed algorithm based on CEEMD is basically consistent with the original spectrum after noise reduction. However, the T2 spectrum curves obtained by EMD and wavelet analysis are significantly different from the original spectrum.
To further verify the effectiveness of the method, random noise was added to the first three groups of data. Then, the method described in this paper was applied to noise reduction of these data. Table 3 shows the SNR comparison processed by the three methods. As shown in Figure 21-23, the T2 spectrum of signal added noise is obtained by inversion. The algorithm proposed in this paper not only improves the signal-to-noise ratio, but also can basically   restore the original T2 spectrum. However, EMD and wavelet analysis can hardly restore the original T2 spectrum.
The above experimental results show that the proposed method is not only suitable for processing high noise signals, but also for processing low noise signals. When dealing with high noise signals, a higher signal-to-noise ratio can be obtained relative to the EMD and wavelet methods, inverting a more accurate structure. When dealing with low   noise signals, the proposed method can better guarantee the accuracy of signal inversion compared to EMD and wavelet methods.

IV. SUMMARY
In this paper, a complete NMR signal model is established, aiming at noise reduction of echo. Simulation results show that orthogonal modulation filter can effectively reduce echo noise in low SNR environment. In view of the noise reduction problem of echo string, wavelet threshold de-noising is generally favored for its excellent de-noising effect. However, for the NMR signal with noise, when it contains small range of effective information, the wavelet threshold de-noising will suppress most of the noise and remove the effective information with a small range at the same time. CEEMD is an improved algorithm of EMD, which not only retains the advantage of EMD in processing non-stationary signals, but also can effectively overcome the modal aliasing problem of EMD. Compared with the heuristic wavelet threshold de-noising method based on deterministic algorithm and the threshold de-noising algorithm based on EMD, the improved CEEMD threshold de-noising method proposed in this paper can not only improve the signal-to-noise ratio, but also get more accurate T2 inversion spectrum while improving SNR.
The proposed CEEMD threshold de-noising process is combined with noise modeling. Fractional Gaussian noise, as the generalization of white noise, reflects the correlation characteristic of noise, well suitable for modeling NMR signals. The standard deviations of noise in the first two IMFs are estimated by the robust estimator. Then, the noise variances in other IMFs can be derived according to the variance relation among the IMFs decomposed from FGn. The estimated variances of noise and the corresponding thresholds are IMF order-dependent. The thresholds are evaluated by the estimated noise variances and the interval threshold scheme is adopted in the proposed CEEMD threshold de-noising. According to the analysis in section 3, the proposed method is not only suitable for processing high noise signals, but also for processing low noise signals. Therefore, this paper provides a strong support for the application of low-field NMR technology in the environment of strong noise. However, this method cannot be used for real-time online systems at this stage. It depends on its complexity. In order to be applied to real-time signal processing, the structure of the algorithm should be further designed in the application. This is also our future research direction.