Method for random noise reduction from MRS oscillating signal using joint algorithms of EMD and TFPF

Abstract

The instant invention relates to a method for noise reduction from a magnetic resonance sounding (MRS) oscillating signal, and more particularly, to a data processing method for reducing random noise contained in MRS oscillating signal based on joint algorithm principles of EMD and TFPF. A MRS oscillating signal is decomposed into different eigen-mode components by using decomposition characteristic of EMD algorithm; then a signal-dominated eigen-mode component is encoded as an instantaneous frequency of an analytical signal of unit amplitude using TFPF algorithm; and random noise is suppressed with the characteristics that the time-frequency distribution of the analytical signal is concentrated along with the instantaneous frequency. The method requires fewer filtering constraints and is simple to operate without need of designing a filtering interval in the time-frequency domain, and has good adaptability to the MRS oscillating signal with a low signal-noise-ratio.

Claims

1. A method for reducing random noise from a magnetic resonance sounding (MRS) oscillating signal for groundwater detection using joint algorithms of Empirical Mode Decomposition (EMD) and Time-Frequency Peak Filtering (TFPF), comprising the following steps: a). performing an EMD decomposition for a MRS observed signal X(n) which is collected by a MRS groundwater detector, to decompose i different eigen-mode components C.sub.1(n), . . . , C.sub.i(n) and one trend term R.sub.i(n) from a high frequency to a low frequency, and X(n)=C.sub.1(n)+. . . +C.sub.i(n)+R.sub.i(n); where n is a discrete sample point; b). performing a Fourier Transform and an auto-correlation analysis for the eigen-mode components C.sub.1(n), . . . , C.sub.i(n) and the trend term R.sub.i(n), synthesizing the Fourier Transform results and the auto-correlation analysis results to obtain noise-dominated mode components and signal-dominated mode components, and extracting the signal-dominated mode components c.sub.1(n), . . . , c.sub.j(n); where j is the number of the signal-dominated mode components; c). using a TFPF algorithm to process the signal-dominated mode components c.sub.1(n), . . . , c.sub.j(n) respectively, to eliminate random noise contained in each signal-dominated mode component to obtain signal components s.sub.1(n), . . . , s.sub.j(n); and d). adding the signal components s.sub.1(n), . . . , s.sub.j(n) to obtain a final MRS oscillating signal s(n) having no random noise interference, s(n)=s.sub.1(n)+. . . +s.sub.j(n), thereby the final MRS oscillating signal with an increased signal-to-noise ratio for groundwater detection is generated.

2. The method of claim 1, wherein the EMD decomposition in the step a) uses EMD algorithm, comprising: first step, identifying all the maximum points and all the minimum points of the observed signal X(n), and performing a cubic spline interpolation for all the maximum points and all the minimum points respectively, to obtain an upper envelope curve E.sub.max(n) and a lower envelope curve E.sub.min(n) of the data; second step, calculating a corresponding average value of the upper and lower envelope curves to obtain a mean curve F.sub.1(n); $F_1\left(n\right)=\frac{E_{\max }\left(n\right)+E_{\min }\left(n\right)}{2}$ third step, subtracting the average value F.sub.1(n) from the observed signal X(n) to obtain a detail component H.sub.1(n), and determining whether the detail component is an eigen-mode function on condition that: {circle around (1)}. a mean of the function is zero; the function is local symmetry, and the number of zero -crossing points and the number of extreme points are the same or differ at most by one; and {circle around (2)}. the sum of the upper and lower envelope values of the function is always zero at any time; if the conditions are not satisfied, repeating the first and second steps for H.sub.1(n) until satisfying the conditions, and obtaining a first eigen-mode component C.sub.1(n); fourth step, subtracting the first eigen-mode component C.sub.1(n) from the observed signal X(n) to obtain an amount of residual R.sub.1(n), that is:
X(n)C.sub.1(n)=R.sub.1(n) repeating the first, second and third steps for the amount of residual R.sub.1(n) as a new signal to be processed, for continuous screening signal; until obtaining a second eigen-mode component C.sub.2(n), and so on, obtaining i different eigen-mode components C.sub.1(n), . . . , C.sub.i(n) and one trend item R.sub.i(n):
R.sub.1(n)C.sub.2(n)=R.sub.2(n), . . . , R.sub.i1(n)C.sub.i(n)=R.sub.i(n)
X(n)=C.sub.1(n)+. . . +C.sub.i(n)+R.sub.i(n).

3. The method of claim 1, wherein using the TFPF algorithm in the step c) comprises: by using the TFPF algorithm, processing the first signal-dominated mode component c.sub.1(n) to obtain a signal component s.sub.1(n) from which the random noise has been eliminated, and the processing comprising: 1). performing an edge data extension for the signal-dominated mode component c.sub.1(n) to obtain an extended signal c.sub.1(m), that is, $c_1^{}\left(m\right)=\{\begin{array}{c}{c}_1\left(1\right),1mp\\ c_1\left(m-p\right),p+1mp+N\\ c_1\left(N\right),p+N+1m2p+N\end{array}$ where N is a length of the original observed signal, p is the number of data points extended at both ends; 2). scaling the extended signal c.sub.1(m): $d_1^{}\left(m\right)=\left(a-b\right).Math.\frac{c_1^{}\left(m\right)-\min \left[c_1^{}\left(m\right)\right]}{\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]}+b$ where d.sub.1(m) is the scaled signal, the coefficients a and b are the maximum value and the minimum value of the scaled signal respectively, satisfying 0.5a=max[d.sub.1(m)], b=min [d.sub.1(m)]0; 3). performing frequency modulation encoding on the scaled signal d.sub.1(m) to obtain an analytical signal of unit amplitude: $z_1\left(m\right)=e^{j2\underset{=0}{\overset{m}{.Math.}}{d_1}^{}\left(\right)}$ where z.sub.1(m) is the analytical signal obtained after the frequency modulation encoding; 4). performing Discrete Pseudo-Wigner-Ville Transform for the obtained analytical signal z.sub.1(m) to calculate a time-frequency distribution of z.sub.1(m): $W_{z_1}\left(m,k\right)=2\underset{l=-L}{\overset{L}{.Math.}}w\left(l\right)z_1\left(m+l\right)z_1^{*}\left(m-l\right)e^{-j4\mathrm{kl}}$ where w(l) is a window function and its width is 2L +1; 5). taking the peak of w.sub.z.sub.1(m,k) as a valid estimate of the signal, that is: .sub.1(m)=arg.sub.kmax [W.sub.z.sub.1(m,k)] where arg.sub.kmax [] represents taking the maximum operator along frequency; 6). inverse-scaling the estimated value of the valid signal, and restoring the signal amplitude: $\left(m\right)=\frac{\left({}^{}\left(m\right)-b\right)\left(\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]\right)}{a-b}+\min \left[c_1^{}\left(m\right)\right]$ 7). removing the extended edge obtained after filtering, to obtain the signal component s.sub.1(n) from which the random noise has been eliminated:
s.sub.1(n)=.sub.1(n+p), 1nN 8). processing the signal-dominated mode components c.sub.2(n), . . . , c.sub.j(n) sequentially in accordance with the above steps 1) to 7), to obtain the signal components s.sub.2(n), . . . , s.sub.j(n) from which the random noise has been eliminated.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flow chart of a method of random noise reduction from a MRS oscillating signals using joint algorithms of EMD and TFPF;

(2) FIG. 2 shows a MRS observed signal containing a random noise;

(3) FIG. 3 shows the eigen-mode components obtained after decomposing the observed signal by using EMD algorithm;

(4) FIG. 4 shows the Fourier Transform of eigen-mode components;

(5) FIG. 5 shows the auto-correlation analysis of eigen-mode components;

(6) FIG. 6 shows the signal-dominated eigen-mode components;

(7) FIG. 7 shows the signal component obtained after reducing the random noise by using TFPF algorithm; and

(8) FIG. 8 shows the MRS oscillating signal obtained after reducing the random noise.

DETAILED DESCRIPTION

(9) The present invention will be now further described in detail below with reference to the accompanying drawings and embodiments.

(10) In accordance with an embodiment of the present invention, a method of random noise reduction from a MRS oscillating signal using joint algorithms of EMD and TFPF, comprises the following steps:

(11) a). performing an EMD decomposition for a MRS observed signal X(n) which is collected by a MRS groundwater detector, to decompose i different eigen-mode components C.sub.1(n), . . . , C.sub.i(n) and one trend term R.sub.i(n) from a high frequency to a low frequency, and X(n)=C.sub.1(n)+ . . . +C.sub.i(n)+R.sub.i(n) ; where n is a discrete sample point;

(12) performing a Fourier Transform and an auto-correlation analysis for the eigen-mode components C.sub.1(n), . . . , C.sub.i(n) and the trend term R.sub.i(n) , synthesizing the Fourier Transform results and the auto-correlation analysis results to obtain noise-dominated mode components and signal-dominated mode components, and extracting the signal-dominated mode components c.sub.1(n), . . . , c.sub.j(n) ; where j is the number of the signal-dominated mode components;

(13) c). using a TFPF algorithm to process the signal-dominated mode components c.sub.1(n), . . . , c.sub.j(n) respectively, to eliminate random noise contained in each signal-dominated mode component to obtain signal components s.sub.1(n), . . . , s.sub.j(n); and

(14) d). adding the signal components s.sub.1(n), . . . , s.sub.j(n) to obtain a final MRS oscillating signal s(n) having no random noise interference, s(n)=s.sub.1(n)+. . . +s.sub.j(n).

(15) The EMD decomposition in the step a) uses EMD algorithm; and the use of the EMD algorithm comprises:

(16) first step, identifying all the maximum points and all the minimum points of the observed signal X(n) , and performing a cubic spline interpolation for all the maximum points and all the minimum points respectively, to obtain an upper envelope curve E.sub.max(n)and a lower envelope curve E.sub.min(n)of the data;

(17) second step, calculating a corresponding average value of the upper and lower envelope curves to obtain a mean curve F.sub.1(n);

(18) $F_1\left(n\right)=\frac{E_{\max }\left(n\right)+E_{\min }\left(n\right)}{2}$

(19) third step, subtracting the average value F.sub.1(n) from the observed signal X(n) to obtain a detail component H.sub.1(n), and determining whether the detail component is an eigen-mode function on condition that:

(20) {circle around (1)}. a mean of the function is zero; the function is local symmetry, and the number of zero-crossing points and the number of extreme points are the same or differ at most by one; and;

(21) {circle around (2)}. the sum of the upper and lower envelope values of the function is always zero at any time;

(22) if the conditions are not satisfied, repeating the first and second steps for H.sub.1(n) until satisfying the conditions, and obtaining a first eigen-mode component C.sub.1(n);

(23) fourth step, subtracting the first eigen-mode component C.sub.1(n) from the observed signal X(n) to obtain an amount of residual R.sub.1(n), that is:
X(n)C.sub.1(n)=R.sub.1(n)

(24) repeating the first, second and third steps for the amount of residual R.sub.1(n) as a new signal to be processed, for continuous screening signal; until obtaining a second eigen-mode component C.sub.2(n), and so on, obtaining i different eigen-mode components C.sub.1(n), C.sub.i(n) and one trend item R.sub.i(n):
R.sub.1(n)C.sub.2(n)=R.sub.2(n) , . . . , R.sub.i1(n)C.sub.i(n)=R.sub.i(n)
X(n)=C.sub.1(n)+. . . +C.sub.i(n)+R.sub.i(n);

(25) The TFPF algorithm in the step c) is used for processing signal-dominated mode components, to obtain the signal components from which the random noise has been eliminated; and the processing comprises:

(26) 1). performing an edge data extension for the signal-dominated mode component c.sub.1(n) to obtain an extended signal c.sub.1(m), that is,

(27) $c_1^{}\left(m\right)=\{\begin{array}{c}{c}_1\left(1\right),1mp\\ c_1\left(m-p\right),p+1mp+N\\ c_1\left(N\right),p+N+1m2p+N\end{array}$
where N is a length of the original observed signal, and p is the number of data points extended at both ends;

(28) 2). scaling the extended signal c.sub.1(m):

(29) $d_1^{}\left(m\right)=\left(a-b\right).Math.\frac{c_1^{}\left(m\right)-\min \left[c_1^{}\left(m\right)\right]}{\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]}+b$
where d.sub.1(m) is the scaled signal, the coefficients a and b are the maximum value and the minimum value of the scaled signal respectively, satisfying 0.5a=max[d.sub.1(m)], b=min[d.sub.1(m)]0;

(30) 3). performing frequency modulation encoding on the scaled signal d.sub.1(m) to obtain an analytical signal of unit amplitude:

(31) $z_1\left(m\right)=e^{j2\underset{=0}{\overset{m}{.Math.}}{d_1}^{}\left(\right)}$
where z.sub.1(m) is the analytical signal obtained after the frequency modulation encoding;

(32) 4). performing Discrete Pseudo-Wigner-Ville Transform for the obtained analytical signal z.sub.1(m) to calculate a time-frequency distribution of z.sub.1(m):

(33) 0$W_{z_1}\left(m,k\right)=2\underset{l=-L}{\overset{L}{.Math.}}w\left(l\right)z_1\left(m+l\right)z_1^{*}\left(m-l\right)e^{-j4\mathrm{kl}}$
where w(l) is a window function and its width is 2L+1;

(34) 5). taking the peak of W.sub.z.sub.1(m,k) as a valid estimate of the signal, that is:
(m)=arg.sub.kmax[W.sub.z.sub.1(m,k)] where arg.sub.kmax[] represents taking the maximum operator along frequency;

(35) 6). inverse-scaling the estimated value of the valid signal, and restoring the signal amplitude:

(36) $\left(m\right)=\frac{\left({}^{}\left(m\right)-b\right)\left(\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]\right)}{a-b}+\min \left[c_1^{}\left(m\right)\right]$

(37) 7). removing the extended edge obtained after filtering, to obtain the signal component s.sub.1(n) from which the random noise has been eliminated:
s.sub.1(n)=(n+p), 1nN

(38) 8). Next, processing the signal-dominated mode components c.sub.2(n), . . . , c.sub.j(n) sequentially in accordance with the above steps 1) to7), to obtain the signal components s.sub.2(n), . . . , S.sub.j(n) from which the random noise has been eliminated.

(39) The embodiment of the method of random noise reduction from a MRS oscillating signal using joint algorithms of EMD and TFPF will be further described in detail:

(40) 1). As shown in FIG. 1, the MRS observed signal X(n) is read; and the observed signal is interfered by a random noise, and is collected by a MRS groundwater detector. FIG. 2 shows the observed signal.

(41) 2). All of extreme points of the observed signal X(n), including the maximum points and the minimum points, are identified; and the cubic spline interpolation is performed for all of the maximum points and the minimum points respectively to obtain the upper envelope curve E.sub.max(n) and the lower envelope curve E.sub.min(n) of the data.

(42) 3). The corresponding average value of the upper and lower envelope curves is calculated to obtain the mean curve F.sub.1(n),

(43) $F_1\left(n\right)=\frac{E_{\max }\left(n\right)+E_{\min }\left(n\right)}{2}.$

(44) 4). The average value F.sub.1(n) is subtracted from the observed signal X(n) to obtain the detail component H.sub.1(n), (H.sub.1(n)=X(n)F.sub.1(n)); it is determined whether the detail component is an eigen-mode function or not, on the following conditions that:

(45) first, the mean of the function is zero; the function is local symmetry; and the number of zero crossing points and the number of extreme points are the same or differ at most by one; and

(46) second, the sum of the upper and lower envelope values is always zero at any time.

(47) If the conditions are not satisfied, steps 2) and 3) are repeated for H.sub.1(n) until the conditions are satisfied, and the first eigen-mode component C.sub.1(n) is obtained.

(48) 5). The first eigen-mode component C.sub.1(n) is subtracted from the observed signal X(n) to obtain the amount of residual R.sub.1(n), i.e., X(n)C.sub.1(n)=R.sub.1(n). The steps 2) 4) are executed repeatedly for the amount of residual R.sub.1(n) as the new signal to be processed for continuous screening signal, until the second eigen-mode component C.sub.2(n) is obtained. And so on, as shown in FIG. 3, fourteen different eigen-mode components C.sub.1(n), . . . , C.sub.14(n) are obtained.

(49) 6). The Fourier transform and the auto-correlation analysis are performed for the eigen-mode components C.sub.1(n), . . . , C.sub.14(n) shown in FIG. 3; and the spectrum of the eigen-mode components can be obtained by using fft command in Matlab, as shown in FIG. 4. The normalized auto-correlation function of the eigen-mode components is calculated by using xcorr command, as shown in FIG. 5. The results of the Fourier transform and the auto-correlation analysis are synthesized to obtain noise-dominated mode components and signal-dominant mode components; and three noise-dominated mode components c.sub.1(n), c.sub.2(n) and c.sub.3(n) are extracted, as shown in FIG. 6.

(50) 7). The first signal-dominated mode component c.sub.1(n) is processed by using the TFPF algorithm; and then the mode component can be encoded as the instantaneous frequency of the analytical signal. Since the edge energy of the time-frequency distribution of the analytical signal is not concentrated, the phenomenon of end-point distortion will be generated when the signal is recovered from the peak of time-frequency distribution. Consequently, it is necessary to perform the edge data extension firstly, to obtain the extended signal c.sub.1(m), that is

(51) $c_1^{}\left(m\right)=\{\begin{array}{c}{c}_1\left(1\right),1mp\\ c_1\left(m-p\right),p+1mp+N\\ c_1\left(N\right),p+N+1m2p+N\end{array}$ where N is the length of the original observed signal, and p is the number of data points extended at both ends.

(52) 8). To avoid a frequency overflow, when the signal is modulated as the frequency information, it is necessary to ensure that the frequency information is within a limited frequency range, consequently the signal c.sub.1(m) needs to be scaled:

(53) $d_1^{}\left(m\right)=\left(a-b\right).Math.\frac{c_1^{}\left(m\right)-\min \left[c_1^{}\left(m\right)\right]}{\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]}+b$ where d.sub.1(m) is the scaled signal, the coefficients a and b are the maximum and minimum values of the scaled signal respectively, satisfying 0.5a=max[d.sub.1(m)], b=min[d.sub.1(m)]0.

(54) 9). The frequency modulation encoding is performed for the scaled signal d.sub.1(m), to obtain the analytical signal of unit amplitude:

(55) $z_1\left(m\right)=e^{j2\underset{=0}{\overset{m}{.Math.}}{d_1}^{}\left(\right)}$
Where z.sub.1(m) is the analytical signal obtained after the frequency modulation encoding.

(56) 10). Discrete Pseudo-Wigner-Ville Transform is performed for the obtained analytical signal z.sub.1(m); and the time-frequency distribution of z.sub.1(m) is calculated:

(57) $W_{z_1}\left(m,k\right)=2\underset{l=-L}{\overset{L}{.Math.}}w\left(l\right)z_1\left(m+l\right)z_1^{*}\left(m-l\right)e^{-j4\mathrm{kl}}$ where w(l) is a window function and its width is 2L+1. After the Wigner-Ville Transform, the signal can be regarded as the distribution of its energy in the joint time domain and frequency domain, which is an important tool for analyzing a non-stationary time variant signal. When the signal is buried in the noise and is a linear function of time, an unbiased estimate of the signal can be obtained. Since the MRS oscillating signal is not a linear function of time, a local linearization processing is performed according to a method of instantaneous frequency estimation, by using a windowed Wigner-Ville Transform, i.e., the Pseudo-Wigner-Ville Transform, to calculate the time-frequency distribution, so that the signal approximately satisfies the condition of the linear instantaneous frequency within a window length; and the deviation caused by non-linearity of the signal can be reduced.

(58) 11). The peak of W.sub.z.sub.1(m, k) is taken as a valid estimate of the signal, that is
(m)=arg.sub.kmax[W.sub.z.sub.1(m,k)]

(59) where arg.sub.kmax[] represents taking the maximum operator along frequency.

(60) 12). The estimated value of the valid signal is inverse-scaled to restore the signal amplitude.

(61) $\left(m\right)=\frac{\left({}^{}\left(m\right)-b\right)\left(\max \left[c_1^{}\left(m\right)\right]-\min \left[c_1^{}\left(m\right)\right]\right)}{a-b}+\min \left[c_1^{}\left(m\right)\right]$

(62) 13). The resulting extended edge obtained after filtering is removed to obtain the signal component s.sub.1(n) from which the random noise has been eliminated.
s.sub.1(n)=(n+p),1nN

(63) 14). Next, the signal-dominated mode components c.sub.2(n) and c.sub.3(n) are sequentially processed in accordance with the above steps 7) to 13), to obtain the signal components from which the random noise has been eliminated, as shown in FIG. 7.

(64) 15). The processed signal components are added together to obtain a final MRS oscillating signal in which the random noise interference is eliminated, as shown in FIG. 8.

(65) Although the invention herein has been described with reference to particular embodiments, it is to be understood that these embodiments are merely illustrative of the principles and applications of the present invention. It is therefore to be understood that numerous modifications may be made to the illustrative embodiments and that other arrangements may be devised without departing from the spirit and scope of the present invention as defined by the following claims.