CN104459809A - Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis - Google Patents

Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis Download PDF

Info

Publication number
CN104459809A
CN104459809A CN201410611932.XA CN201410611932A CN104459809A CN 104459809 A CN104459809 A CN 104459809A CN 201410611932 A CN201410611932 A CN 201410611932A CN 104459809 A CN104459809 A CN 104459809A
Authority
CN
China
Prior art keywords
signal
mrs
independent component
component analysis
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201410611932.XA
Other languages
Chinese (zh)
Other versions
CN104459809B (en
Inventor
田宝凤
周媛媛
王悦
李振宇
谢燕妮
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN201410611932.XA priority Critical patent/CN104459809B/en
Publication of CN104459809A publication Critical patent/CN104459809A/en
Application granted granted Critical
Publication of CN104459809B publication Critical patent/CN104459809B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

The invention discloses a full-wave nuclear magnetic resonance signal denoising method based on independent component analysis. The full-wave nuclear magnetic resonance signal denoising method is mainly used for power frequency harmonic interference or certain single-frequency interference in full-wave nuclear magnetic resonance signals. The method comprises the steps that firstly, utilizing a nuclear magnetic resonance sounding water detecting instrument for collecting MRS signals, obtaining the frequency of power frequency harmonic interference or certain single-frequency interference contained in the collected signals through spectral analysis, and adopting a digital orthogonal method for constructing input channel signals to solve the underdetermined blind source separation problem; then, using the constructed input channel signals and the collected MRS signals as input signals for conducting independent component analysis to obtain separated MRS signals; finally, adopting the spectrum correcting method for solving the problem that the separated MRS signals have uncertain amplitude, and then extracting denoised MRS signals. According to the full-wave nuclear magnetic resonance signal denoising method, the problem of the underdetermined blind source and the problem of the uncertain amplitude usually existing in independent component analysis are solved, the power frequency harmonic interference or the certain single-frequency interference in the noised MRS signals is effectively filtered out, and compared with a traditional MRS signal denoising method, the full-wave nuclear magnetic resonance signal denoising method has the advantages of being high in arithmetic speed, high in signal to noise ratio, high in practicability and the like.

Description

A kind of all-wave NMR signal noise filtering method based on independent component analysis
Technical field
The present invention relates to a kind of nuclear magnetic resonance depth measurement (Magnetic Resonance Sounding, MRS) signal noise filtering method, specifically utilize the filtering method of the all-wave NMR signal noise based on independent component analysis principle.
Background technology
Nuclear magnetic resonance depth measurement (Magnetic Resonance Sounding, MRS) underground water direct detection method that is a kind of high-new, harmless, high resolution that technology is that developed recently gets up, it is a kind of Weak Signal Detection, and its ultimate principle is that the NMR signal produced by Hydrogen Proton resonant transition in Underground water realizes underground water detection.The general expression of all-wave NMR signal is:
(10)
Four key parameter: E of MRS signal are comprised in formula (10) 0(q), ω 0, represent the angular Larmor frequency from the signal initial amplitude that stratum water cut is correlated with, the relaxation time (the relaxation time scope of free water be 30ms-1000ms) relevant with water-bearing zone factor of porosity, zones of different with different numerical value, the initial phase relevant with water-bearing zone electric conductivity respectively.
Although MRS method has all advantages in underground water detection, detect nuclear magnetic signal at present and belong to and receive volt level signal, extremely faint, be very easily subject to the impact of neighbourhood noise.Especially when the Hz noise in neighbourhood noise or a certain mono-tone interference frequency are near Larmor frequency, the complete None-identified of MRS signal will be caused, thus be difficult to accurately extract the key parameter such as initial amplitude, relaxation time, have a strong impact on follow-up inversion interpretation work.
At present, for the Theories and methods of NMR signal noise filtering, domestic and international experts and scholars have carried out a large amount of research work.On hardware implementing, generally adopt Trushkin D V, Shushakov O A, the data investigation method that Legchenko A V proposes in " The potential of a noise-reducing antenna for surface NMRgroundwater surveys in the earth's magnetic field " (" Geophysical Prospecting " the 42nd volume 8 phase in 1994: 855-862 page) and Konstantinos Chalikakis, Mette Ryom Nielsen, Anatoly Legchenko is at " MRS applicability for a study of glacialsedimentary aquifers in Central Jutland, Denmark " (" Journal of AppliedGeophysics " the 66th volume 3 phase in 2008: 176 – 187 pages) middle " 8 " font coil raising signal to noise ratio (S/N ratio) proposed, but can cause the reduction of instrument investigation depth and work efficiency low, the hyperchannel GeoMRI instrument of in recent years representative VistaClara company of U.S. development and French IRIS company release last word NUMISPOLY instrument and all propose to adopt reference channel mode, improve signal to noise ratio (S/N ratio) based on adaptive noise cancellation method.The de-noising thought of above-mentioned algorithm is all realize SNR estimation and compensation based on the correlativity of noise in reference channel noise and main channel signals and associated noises, but due to the irregularities of some noise and the uncertainty of instability and mixed mechanism, limit the application of algorithm.
Due in most cases, source signal, noise and mixed mechanism are all unknown, and therefore, a kind of independent component analysis (Independent Component Analysis, ICA) method of rising in recent years solves problems.ICA is the effective separation means of one developed on the basic model of blind source separating problem, and its basic model is X=AS (X is observation signal, and A is unknown linear combination matrix, and S is source signal) and solving model is Y=W t(Y is separation signal to X, A -1=W tfor separation matrix).ICA is transformed into solving model by certain optimized algorithm and objective function from basic model, and then obtains separation signal.
Patent CN101908138A discloses " a kind of identification method of image target of synthetic aperture radar based on noise independent component analysis, belongs to the technical field of diameter radar image target identification "; Patent CN102697493A discloses " a kind of method that in EEG signals fast, eye electricity artefact automatically identifies and removes ", belongs to medical signals field; Patent CN102867189A discloses one " the ADAPTIVE MIXED separation of images method based on independent component analysis ", belongs to image processing field; Patent CN103854660A discloses " a kind of four Mike's sound enhancement methods based on independent component analysis ", belongs to field of voice signal.Visible independent component analysis has been successfully applied the every field of signal transacting, but there is not yet it and be applied in the noise filtering of MRS signal.
At present, when applying traditional independent component analysis and carrying out Signal separator, main Problems existing is: owe to determine blind source separating problem, separation signal amplitude is indefinite.So-called owing determines blind source separating problem, refers in the noise source that MRS signal comprises, and noise is complicated, and cause independent source number unknown, and consider the restriction of instrument and equipment, the number of the observation signal collected is less than the situation of the number of source signal.Problems is undoubtedly by the application of restriction ICA in the filtering of nuclear magnetic resonance depth measurement signal noise.In addition, if noisy MRS signal is after ICA, the amplitude being separated MRS signal can not get accurate estimation, then stratum water cut estimated error rate in follow-up inversion interpretation can be caused to increase, affect engineering progress.
Summary of the invention
Technical matters to be solved by this invention is the filtering method providing a kind of all-wave NMR signal noise based on independent component analysis principle.
The present invention is achieved in that a kind of all-wave NMR signal noise filtering method based on independent component analysis, comprises the following steps:
Step (1): utilize nuclear magnetic resonance depth measurement to visit water instrument and collect one group of observation MRS signal X 1(t);
Step (2): to the observation MRS signal X gathered 1t () carries out Fourier transform, obtain its frequency spectrum, determines the frequency f of industrial frequency harmonic or a certain mono-tone interference N (t) contained in signal 0, f 1f n;
What industrial frequency harmonic interference represented is the integral multiple of 50Hz, and such as can represent industrial frequency harmonic a: X=a*cos (2*pi*f*t) as follows, a is wherein amplitude, can get arbitrary constant, and f is 2350Hz (being 47 times of 50); A certain mono-tone interference refers to the interference that frequency is not the integral multiple of 50Hz.
Step (3): according to the industrial frequency harmonic determined in step (2) or a certain mono-tone interference N (t) frequency f 0, f 1f n, adopt digital quadrature method structure input channel signal X 2(t), X 3(t) ... X 2n(t), X 2n+1(t);
Step (4): by the input channel signal X of structure in step (3) 2(t), X 3(t) ... X 2n(t), X 2n+1(t) and observation MRS signal X 1t () is in the lump as input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)] carry out independent component analysis, obtain separation signal
Step (5): separation signal is carried out spectrum analysis respectively, according to desirable MRS signal S (t) feature determination target MRS signal
Step (6): adopt spectrum correction method to the separation signal after independent component analysis carry out spectrum recovery, obtain MRS signal Y (t) after denoising.
Further, the concrete steps of the digital quadrature method in described step (3) are:
First, observation signal X is supposed 1t () is made up of, represented by (1): X desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t) 1(t)=AS (t)+BN (t)=AS (t)+Bsin (2 π f 0t+ θ 0), wherein A, B are respectively two arbitrary constants, f 0for the industrial frequency harmonic that obtains through Frequency Estimation method or a certain mono-tone interference frequency, θ 0for the phase place of industrial frequency harmonic interference;
Secondly, due to phase theta 0be difficult to estimate, therefore formula (1) done further decomposition:
X 1(t)=A·S(t)+B·sin(2πf 0t+θ 0)=
(2)
A·S(t)+B·sin(2πf 0t)·cos(θ 0)+B·cos(2πf 0t)·sin(θ 0)
Cos (θ in formula (2) 0), sin (θ 0) be time-independent constant, if:
B 1=B·cos(θ 0),B 2=B·sin(θ 0)
(3)
Then:
X 1(t)=A·S(t)+B 1·sin(2πf 0t)+B 2·cos(2πf 0t) (4)
Finally, by observation signal X 1t () represents the form of an accepted way of doing sth (4), therefore construct input channel signal:
X 2(t)=B 1sin (2 π f 0t), X 3(t)=B 2cos (2 π f 0t) (5) are as observation signal X 1t () be two-way input channel signal in addition, if when having the interference of multiple industrial frequency harmonic or a certain mono-tone interference, can carry out according to the method described above constructing multiple input channel signal equally, that is: X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)];
Further, the independent component analysis in described step (4), adopts the fast fixed point FastICA algorithm based on negentropy to carry out noisy MRS signal X 1the separation of (t), concrete steps are:
1), digital quadrature method is adopted to obtain input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)], carried out pre-service, comprise decentralization and albefaction, object makes the average of observation signal be zero and correlativity is less, obtains pretreated signal Z (t)=[Z 1(t); Z 2(t); Z 3(t); ... Z 2n(t); Z 2n+1(t)];
2), judge the Gaussian of desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t), be convenient to select suitable nonlinear function to carry out independent component separation, the measure formulas according to kurtosis:
Kurt=E{x 4}/E 2{ x 2} -3(6) kurtosis value that can be calculated two kinds of signals is all less than zero, known S (t) and N (t) all belong to sub-Gaussian signals, therefore, selecting type (7) substitutes into iterative formula (8) as the nonlinear function of FastICA algorithm:
G ( Y ) = 1 4 Y 4 - - - ( 7 )
g(Y)=Y 3
Wherein, g () the single order derived function that is G (); Y=(w i t(p) X (t));
3), the condition of convergence is set, and initialization weight vector;
4), according to negentropy maximization condition, the iterative formula (8) adopting Newton iteration method to obtain and normalization formula (9) solve separation matrix W, thus obtain separation signal
w i(p+1)=E{Xg(w i T(p)X(t))}-E{G″(w i T(p)X(t))}w i(p) (8)
w i(p+1)=w i(p+1)/||w i(p+1)|| (9)
P (p=1,2,3...n) in formula (9) is iterative steps, w ip () is one of them row vector corresponding to W; Formula (9) realizes w i(p+1) normalization.
Further, the spectrum correcting method in described step (6), concrete steps are:
1), independent component analysis is separated the target MRS signal obtained carry out Fourier transform with desirable MRS signal S (t), obtain both frequency spectrums;
2), contrast both the difference degree of spectral magnitude at Larmor frequency place, the former frequency spectrum multiple is returned to the amplitude of the latter's frequency spectrum;
3), by the frequency spectrum of separation signal after correcting through Fourier inversion to time domain, obtain MRS signal Y (t) after denoising.
The present invention compared with prior art, beneficial effect is: the present invention proposes the all-wave NMR signal noise filtering method based on independent component analysis, not only achieve the industrial frequency harmonic interference near MRS signal Larmor frequency or effective filtering of mono-tone interference, and the method utilizes the input channel signal of digital quadrature method structure except observation signal passage, without the need to the requirement of traditionally ICA, extra instrument and equipment is used to gather input multiplexer signal, break the harsh conditions of ICA application, save a large amount of manpower and materials.Meanwhile, the separation signal amplitude indeterminate problem utilizing spectrum correction method to solve traditional IC A to be used to have, for the initial amplitude extracting MRS signal has accurately and efficiently carried out sufficient preparation.Compared with traditional MRS signal antinoise method, the inventive method can reduce MRS signal message and lose, and fast operation, signal to noise ratio (S/N ratio) is high.In addition, the inventive method also opens the new world of independent component analysis in NMR signal de-noising field, is the applications expanding thinking of follow-up independent component analysis.
Accompanying drawing explanation
Fig. 1 is based on the FB(flow block) of the all-wave NMR signal noise filtering method of independent component analysis;
Fig. 2 Fixed point algorithm FB(flow block);
Fig. 3 Spectrum Correction process schematic;
The desirable MRS signal of Fig. 4 (a);
The desirable MRS signal spectrum of Fig. 4 (b);
The noisy MRS signal of Fig. 5 (a);
The noisy MRS signal spectrum of Fig. 5 (b);
Fig. 6 (a) is separated MRS signal;
Fig. 6 (b) is separated MRS signal spectrum;
MRS signal after Fig. 7 (a) Spectrum Correction;
MRS signal spectrum after Fig. 7 (b) Spectrum Correction;
Fig. 8 (a) observes MRS signal;
Fig. 8 (b) observes MRS signal spectrum;
Fig. 9 (a) observes MRS separation signal;
Fig. 9 (b) observes MRS separation signal frequency spectrum;
Actual measurement MRS signal before and after Figure 10 (a) ICA process;
Actual measurement MRS signal spectrum before and after Figure 10 (b) ICA process.
Embodiment
In order to make object of the present invention, technical scheme and advantage clearly understand, below in conjunction with drawings and Examples, the present invention is further elaborated.Should be appreciated that specific embodiment described herein only in order to explain the present invention, be not intended to limit the present invention.
As shown in Figure 1, a kind of all-wave NMR signal noise filtering method based on independent component analysis, comprises the following steps:
Step (1): utilize nuclear magnetic resonance depth measurement to visit water instrument and collect one group of observation MRS signal X 1(t);
Step (2): to the observation MRS signal X gathered 1t () carries out Fourier transform, obtain its frequency spectrum, determines the frequency f of industrial frequency harmonic or a certain mono-tone interference N (t) contained in signal 0, f 1f n;
Step (3): according to the industrial frequency harmonic determined in step (2) or a certain mono-tone interference N (t) frequency f 0, f 1f n, adopt digital quadrature method structure input channel signal X 2(t), X 3(t) ... X 2n(t), X 2n+1(t);
Step (4): by the input channel signal X of structure in step (3) 2(t), X 3(t) ... X 2n(t), X 2n+1(t) and observation MRS signal X 1t () is in the lump as input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)] carry out independent component analysis, obtain separation signal
Step (5): separation signal is carried out spectrum analysis respectively, according to desirable MRS signal S (t) feature determination target MRS signal
Step (6): adopt spectrum correction method to the separation signal after independent component analysis carry out spectrum recovery, obtain MRS signal Y (t) after denoising.
A kind of concrete steps of digital quadrature method of the all-wave NMR signal noise filtering method based on independent component analysis are:
First, observation signal X is supposed 1t () is made up of, represented by (1): X desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t) 1(t)=AS (t)+BN (t)=AS (t)+Bsin (2 π f 0t+ θ 0), wherein A, B are respectively two arbitrary constants, f 0for the industrial frequency harmonic that obtains through Frequency Estimation method or a certain mono-tone interference frequency, θ 0for the phase place of industrial frequency harmonic interference,
Secondly, due to phase theta 0be difficult to estimate, therefore formula (1) done further decomposition:
X 1(t)=A·S(t)+B·sin(2πf 0t+θ 0)=
(2)
A·S(t)+B·sin(2πf 0t)·cos(θ 0)+B·cos(2πf 0t)·sin(θ 0)
Cos (θ in formula (2) 0), sin (θ 0) be time-independent constant, if:
B 1=B·cos(θ 0),B 2=B·sin(θ 0)
(3)
Then:
X 1(t)=A·S(t)+B 1·sin(2πf 0t)+B 2·cos(2πf 0t) (4)
Finally, by observation signal X 1t () represents the form of an accepted way of doing sth (4), therefore construct input channel signal:
X 2(t)=B 1sin (2 π f 0t), X 3(t)=B 2cos (2 π f 0t) (5) are as observation signal X 1t () be two-way input channel signal in addition, if when having the interference of multiple industrial frequency harmonic or a certain mono-tone interference, can carry out according to the method described above constructing multiple input channel signal equally, that is: X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)];
As shown in Figure 2, independent composition analysis algorithm provided by the invention, the Fixed point algorithm based on negentropy (FastICA) adopting the people such as Hyvarinen to propose carries out noisy MRS signal X 1the separation of (t), concrete steps are:
1), digital quadrature method is adopted to obtain input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)], carried out pre-service, comprise decentralization and albefaction, object makes the average of observation signal be zero and correlativity is less, obtains pretreated signal Z (t)=[Z 1(t); Z 2(t); Z 3(t); ... Z 2n(t); Z 2n+1(t)];
2), judge the Gaussian of desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t), be convenient to select suitable nonlinear function to carry out independent component separation, the measure formulas according to kurtosis:
Kurt=E{x 4}/E 2{ x 2} -3(6) kurtosis value that can be calculated two kinds of signals is all less than zero, known S (t) and N (t) all belong to sub-Gaussian signals, therefore, selecting type (7) substitutes into iterative formula (8) as the nonlinear function of FastICA algorithm:
G ( Y ) = 1 4 Y 4 - - - ( 7 )
g(Y)=Y 3
Wherein, g () the single order derived function that is G (); Y=(w i t(p) X (t));
3), the condition of convergence is set, and initialization weight vector;
4), according to negentropy maximization condition, the iterative formula (8) adopting Newton iteration method to obtain and normalization formula (9) solve separation matrix W, thus obtain separation signal
w i(p+1)=E{Xg(w i T(p)X(t))}-E{G″(w i T(p)X(t))}w i(p) (8)
W i(p+1)=w i(p+1)/|| w i(p+1) || the p (p=1,2,3...n) in (9) formula (9) is iterative steps, w ip () is one of them row vector corresponding to W;
Formula (9) realizes w i(p+1) normalization, wherein needs in separation matrix W to judge | w i(p)-w i(p+1) | whether being less than the condition of convergence, in this way, then stopping iteration, if not being return iterative formula (8), re-starting orthogonal, normalization.
As shown in Figure 3, spectrum correcting method provided by the invention, the spectrum correcting method in described step (6), concrete steps are:
1), independent component analysis is separated the target MRS signal (signal as shown in Fig. 3 B1) obtained and desirable MRS signal S (t) (signal as shown in Fig. 3 A1) carries out Fourier transform, obtain both frequency spectrums, signal collection of illustrative plates as shown in Fig. 3 B2 is the frequency spectrum of target MRS signal, the frequency spectrum of desirable MRS signal as shown in figure 3 a 2;
2), contrast both the difference degree of spectral magnitude at Larmor frequency place, the former frequency spectrum multiple is returned to the amplitude of the latter's frequency spectrum;
3), by the frequency spectrum (frequency spectrum of separation signal as shown in Figure 3 C) of separation signal after correcting through Fourier inversion to time domain, obtain MRS signal Y (t) after denoising (the MRS signal after denoising as shown in Figure 3 D).
Embodiment 1
The present embodiment is the emulation experiment of the inventive method of carrying out under MATLAB 7.0 programmed environment.
Based on the simulation algorithm of the all-wave NMR signal noise filtering method of independent component analysis, with reference to Fig. 1, comprise the following steps:
Step (1): utilizing formula (10) to construct Larmor frequency is 2325Hz, amplitude e 0for 150nV, the relaxation time for the desirable MRS signal of 0.15s, as shown in Fig. 4 (a), the frequency spectrum that Fig. 4 (b) is its correspondence.Near this signal Larmor frequency, add the Hz noise of 2350Hz and 2300Hz, form through certain linear combination the observation MRS signal X that signal to noise ratio (S/N ratio) is-6.7471dB 1t () (being row vector), as shown in Fig. 5 (a);
Step (2): to observation MRS signal X 1t () carries out Fourier transform, obtain its frequency spectrum, as shown in Fig. 5 (b).Determine that the work frequency in signals and associated noises is respectively f 0=2300Hz, f 1=2350Hz;
Step (3): according to the noise frequency determined in step (2), adopts digital quadrature method structure input channel signal X 2(t)=sin (2 π f 0t), X 3(t)=cos (2 π f 0t), X 4(t)=sin (2 π f 1t), X 5(t)=cos (2 π f 1t);
Step (4): by the input channel signal X of structure in step (3) 2(t), X 3(t), X 4(t), X 5(t) and observation MRS signal X 1t () is in the lump as input signal X (t)=[X 1(t); X 2(t); X 3(t); X 4(t); X 5(t)] carry out independent component analysis, obtain independent component and separation signal
Step (5): by separation signal carry out spectrum analysis respectively, according to desirable MRS signal S (t) feature determination target MRS signal as shown in Fig. 6 (a), the frequency spectrum that Fig. 6 (b) is its correspondence, its Time Domain Amplitude is about 2.5nV as seen, has a long way to go with pure MRS signal amplitude;
Step (6): adopt spectrum correction method to the separation signal after independent component analysis carry out spectrum recovery, obtain MRS signal Y (t) after denoising, as shown in Fig. 7 (a), the frequency spectrum that Fig. 7 (b) is its correspondence;
In order to verify the practicality of the inventive method, MRS signal Y (t) after denoising being carried out signal to noise ratio (S/N ratio) (SNR) and having estimated.As calculated, its SNR=23.2521dB, the SNR before being comparatively separated improves 29.9992 dB; Then envelope extraction and data fitting have been carried out to Y (t), to obtain the key parameter initial amplitude E of separation signal 0(q) and relaxation time T 2 *, can be calculated, E 0(q)=152.3249nV, T 2 *=0.1453s, relative error is respectively 1.55% ,-3.11%, all controls within ± 5%, meets application requirement.
Embodiment 2
The MRS signal that the present embodiment gathers using enamelware pot town, Changchun (this ground Larmor frequency is for 2326Hz) is as the handling object of the inventive method.As shown in Figure 1, based on the all-wave NMR signal noise filtering method of independent component analysis, comprise the following steps:
Step (1): utilize nuclear magnetic resonance depth measurement (MRS) to visit water instrument and collect one group of observation MRS signal X 1(t) (being row vector), as shown in Fig. 8 (a), calculating its signal to noise ratio (S/N ratio) is SNR=-4.9186dB;
Step (2): to the observation MRS signal X gathered 1t () carries out Fourier transform, obtain its frequency spectrum, as shown in Fig. 8 (b), can see that this signal is at f 0=2300Hz, f 1=2350Hz, f 2all there is stronger industrial frequency harmonic interference at=2450Hz place, in addition, at f 3=2383.5Hz place receives strong mono-tone interference;
Step (3): according to the noise frequency determined in step (2), adopts digital quadrature method structure input channel signal X 2(t)=sin (2 π f 0t), X 3(t)=cos (2 π f 0t), X 4(t)=sin (2 π f 1t), X 5(t)=cos (2 π f 1t), X 6(t)=sin (2 π f 2t), X 7(t)=cos (2 π f 2t), X 8(t)=sin (2 π f 3t), X 9(t)=cos (2 π f 3t);
Step (4): by the input channel signal X of structure in step (3) 2(t), X 3(t), X 4(t), X 5(t), X 6(t), X 7(t), X 8(t), X 9(t) and observation MRS signal X 1t () is in the lump as input signal X (t)=[X 1(t); X 2(t); X 3(t); X 4(t); X 5(t); X 6(t); X 7(t); X 8(t); X 9(t)] carry out independent component analysis, obtain independent component and separation signal
Step (5): by separation signal carry out spectrum analysis respectively, according to desirable MRS signal S (t) feature determination target MRS signal as shown in Fig. 9 (a), be observation MRS separation signal, the frequency spectrum that Fig. 9 (b) is its correspondence;
Step (6): adopt spectrum correction method to the separation signal after independent component analysis carry out spectrum recovery, obtain MRS signal Y (t) that the signal to noise ratio (S/N ratio) after denoising is 2.9698dB, light-colored part as larger in area in Figure 10 (a) is depicted as raw observation MRS signal, the frequency spectrum that the dark parts that area is less in Figure 10 (a) is the signal to noise ratio (S/N ratio) after denoising be MRS signal Y (t) of 2.9698dB, Figure 10 (b) is its correspondence;
The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, all any amendments done within the spirit and principles in the present invention, equivalent replacement and improvement etc., all should be included within protection scope of the present invention.

Claims (4)

1., based on an all-wave NMR signal noise filtering method for independent component analysis, it is characterized in that, comprise the following steps:
Step (1): utilize nuclear magnetic resonance depth measurement to visit water instrument and collect one group of observation MRS signal X 1(t);
Step (2): to the observation MRS signal X gathered 1t () carries out Fourier transform, obtain its frequency spectrum, determines the frequency f of industrial frequency harmonic or a certain mono-tone interference N (t) contained in signal 0, f 1f n;
Step (3): according to the industrial frequency harmonic determined in step (2) or a certain mono-tone interference N (t) frequency f 0, f 1f n, adopt digital quadrature method structure input channel signal X 2(t), X 3(t) ... X 2n(t), X 2n+1(t);
Step (4): by the input channel signal X of structure in step (3) 2(t), X 3(t) ... X 2n(t), X 2n+1(t) and observation MRS signal X 1t () is in the lump as input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)] carry out independent component analysis, obtain separation signal
Step (5): separation signal is carried out spectrum analysis respectively, according to desirable MRS signal S (t) feature determination target MRS signal
Step (6): adopt spectrum correction method to the separation signal after independent component analysis carry out spectrum recovery, obtain MRS signal Y (t) after denoising.
2., according to a kind of all-wave NMR signal noise filtering method based on independent component analysis according to claim 1, it is characterized in that, the concrete steps of the digital quadrature method in described step (3) are:
First, observation signal X is supposed 1t () is made up of, represented by (1): X desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t) 1(t)=AS (t)+BN (t)=AS (t)+Bsin (2 π f 0t+ θ 0), wherein A, B are respectively two arbitrary constants, f 0for the industrial frequency harmonic that obtains through Frequency Estimation method or a certain mono-tone interference frequency, θ 0for the phase place of industrial frequency harmonic interference;
Secondly, due to phase theta 0be difficult to estimate, therefore formula (1) done further decomposition:
X 1(t)=A·S(t)+B·sin(2πf 0t+θ 0)= (2)
A·S(t)+B·sin(2πf 0t)·cos(θ 0)+B·cos(2πf 0t)·sin(θ 0)
Cos (θ in formula (2) 0), sin (θ 0) be time-independent constant, if:
B 1=B·cos(θ 0),B 2=B·sin(θ 0)
(3)
Then:
X 1(t)=A·S(t)+B 1·sin(2πf 0t)+B 2·cos(2πf 0t) (4)
Finally, by observation signal X 1t () represents the form of an accepted way of doing sth (4), therefore construct input channel signal:
X 2(t)=B 1·sin(2πf 0t),X 3(t)=B 2·cos(2πf 0t) (5)
As observation signal X 1t () be two-way input channel signal in addition, if when having the interference of multiple industrial frequency harmonic or a certain mono-tone interference, can carry out according to the method described above constructing multiple input channel signal equally, that is: X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)].
3. according to a kind of all-wave NMR signal noise filtering method based on independent component analysis according to claim 1, it is characterized in that, independent component analysis in described step (4), adopts the fast fixed point FastICA algorithm based on negentropy to carry out noisy MRS signal X 1the separation of (t), concrete steps are:
1), digital quadrature method is adopted to obtain input signal X (t)=[X 1(t); X 2(t); X 3(t); ... X 2n(t); X 2n+1(t)], carried out pre-service, comprise decentralization and albefaction, object makes the average of observation signal be zero and correlativity is less, obtains pretreated signal Z (t)=[Z 1(t); Z 2(t); Z 3(t); ... Z 2n(t); Z 2n+1(t)];
2), judge the Gaussian of desirable MRS signal S (t) and industrial frequency harmonic or a certain mono-tone interference N (t), be convenient to select suitable nonlinear function to carry out independent component separation, the measure formulas according to kurtosis:
kurt=E{x 4}/E 2{x 2} -3(6)
The kurtosis value that can be calculated two kinds of signals is all less than zero, and known S (t) and N (t) all belong to sub-Gaussian signals, and therefore, selecting type (7) substitutes into iterative formula (8) as the nonlinear function of FastICA algorithm:
G ( Y ) = 1 4 Y 4 - - - ( 7 )
g(Y)=Y 3
Wherein, g () the single order derived function that is G (); Y=(w i t(p) X (t));
3), the condition of convergence is set, and initialization weight vector;
4), according to negentropy maximization condition, the iterative formula (8) adopting Newton iteration method to obtain and normalization formula (9) solve separation matrix W, thus obtain separation signal
w i(p+1)=E{Xg(w i T(p)X(t))}-E{G″(w i T(p)X(t))}w i(p) (8)
w i(p+1)=w i(p+1)/||w i(p+1)|| (9)
P (p=1,2,3...n) in formula (9) is iterative steps, w ip () is one of them row vector corresponding to W; Formula (9) realizes w i(p+1) normalization.
4., according to a kind of all-wave NMR signal noise filtering method based on independent component analysis according to claim 1, it is characterized in that, the spectrum correcting method in described step (6), concrete steps are:
1), independent component analysis is separated the target MRS signal obtained carry out Fourier transform with desirable MRS signal S (t), obtain both frequency spectrums;
2), contrast both the difference degree of spectral magnitude at Larmor frequency place, the former frequency spectrum multiple is returned to the amplitude of the latter's frequency spectrum;
3), by the frequency spectrum of separation signal after correcting through Fourier inversion to time domain, obtain MRS signal Y (t) after denoising.
CN201410611932.XA 2014-10-30 2014-10-30 Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis Expired - Fee Related CN104459809B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410611932.XA CN104459809B (en) 2014-10-30 2014-10-30 Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410611932.XA CN104459809B (en) 2014-10-30 2014-10-30 Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis

Publications (2)

Publication Number Publication Date
CN104459809A true CN104459809A (en) 2015-03-25
CN104459809B CN104459809B (en) 2017-04-12

Family

ID=52906153

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410611932.XA Expired - Fee Related CN104459809B (en) 2014-10-30 2014-10-30 Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis

Country Status (1)

Country Link
CN (1) CN104459809B (en)

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777442A (en) * 2015-04-07 2015-07-15 吉林大学 MRS (magnetic resonance sounding) FID (frequency identity) signal noise inhibition method
CN105115622A (en) * 2015-08-12 2015-12-02 合肥工业大学 Denoising algorithm of fiber Raman temperature sensing system based on independent component analysis
CN105549097A (en) * 2015-12-22 2016-05-04 吉林大学 Transient electromagnetic signal power frequency and harmonic interference elimination method and apparatus thereof
CN105954170A (en) * 2016-05-04 2016-09-21 中国石油大学(华东) Nuclear magnetic porosity calculating method considering background noise signal quantity
CN106361328A (en) * 2016-10-21 2017-02-01 电子科技大学 Method for electroencephalogram signal extraction under magnetic resonance environment
CN106931905A (en) * 2017-03-09 2017-07-07 北京理工大学 A kind of digital Moiré patterns phase extraction method based on nonlinear optimization
CN106955112A (en) * 2017-03-17 2017-07-18 泉州装备制造研究所 Brain wave Emotion recognition method based on Quantum wavelet neural networks model
CN107092925A (en) * 2017-03-30 2017-08-25 中国人民解放军国防科学技术大学 Cerebral function magnetic resonance imaging blind source separation method based on SIM algorithms in groups
CN107121705A (en) * 2017-04-28 2017-09-01 中南大学 A kind of ground penetrating radar echo signals Denoising Algorithm compared based on automatic anti-phase correction and kurtosis value
CN108318764A (en) * 2018-03-28 2018-07-24 国网上海市电力公司 A kind of earthing or grounding means shock response test jamproof system and method
CN108846407A (en) * 2018-04-20 2018-11-20 太原理工大学 The nuclear magnetic resonance image classification method of brain network is not known based on independent element high order
CN109143389A (en) * 2018-08-01 2019-01-04 吉林大学 A kind of three-dimensional industrial frequency interference source for nuclear-magnetism quantifies orienting device and measurement method
CN109828318A (en) * 2019-01-25 2019-05-31 吉林大学 A kind of magnetic resonance depth measurement signal noise filtering method based on variation mode decomposition
CN110058312A (en) * 2018-10-22 2019-07-26 南方科技大学 A kind of method, apparatus and terminal device inhibiting the interference of earth magnetism near field noise
CN111368778A (en) * 2020-03-13 2020-07-03 哈尔滨理工大学 Weak signal noise stripping method based on intelligent optimization algorithm
CN112036453A (en) * 2020-08-14 2020-12-04 哈尔滨工程大学 Blind source separation method based on quantum rhinoceros search mechanism
CN112180454A (en) * 2020-10-29 2021-01-05 吉林大学 Magnetic resonance underground water detection random noise suppression method based on LDMM
CN112370066A (en) * 2020-09-30 2021-02-19 北京工业大学 Brain-computer interface method of stroke rehabilitation system based on generation of countermeasure network

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090174402A1 (en) * 2008-01-07 2009-07-09 Baker Hughes Incorporated Joint Compression of Multiple Echo Trains Using Principal Component Analysis and Independent Component Analysis
CN101788656A (en) * 2010-01-29 2010-07-28 东南大学 Method for recognizing function response signal under function nuclear magnetic resonance scan
CN101820338A (en) * 2010-03-30 2010-09-01 中国科学院武汉物理与数学研究所 Method for synchronizing decimation filter of digital receiver of nuclear magnetic resonance spectrometer

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090174402A1 (en) * 2008-01-07 2009-07-09 Baker Hughes Incorporated Joint Compression of Multiple Echo Trains Using Principal Component Analysis and Independent Component Analysis
CN101788656A (en) * 2010-01-29 2010-07-28 东南大学 Method for recognizing function response signal under function nuclear magnetic resonance scan
CN101820338A (en) * 2010-03-30 2010-09-01 中国科学院武汉物理与数学研究所 Method for synchronizing decimation filter of digital receiver of nuclear magnetic resonance spectrometer

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
CHUANDONG JIANG 等: ""Statistical stacking and adaptive notch filter to remove high-level electromagnetic noise from MRS measurements"", 《NEAR SURFACE GEOPHYSICS》 *
YEN-CHIEH OUYANG 等: ""Independent Component Analysis for Magnetic Resonance Image Analysis"", 《EURASIP JOURNAL ON ADVANCES IN SIGNAL PROCESSING》 *
张红娟 等: ""核独立成分分析在fMRI数据中的应用"", 《控制工程》 *
田宝凤 等: ""核磁共振信号工频谐波的自适应滤除方法"", 《吉林大学学报(信息科学版)》 *
魏巍等: ""基于独立分量分析的工频干扰消除技术"", 《计算机应用研究》 *

Cited By (27)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777442B (en) * 2015-04-07 2017-06-16 吉林大学 A kind of nuclear magnetic resonance depth measurement FID signal noise suppressing method
CN104777442A (en) * 2015-04-07 2015-07-15 吉林大学 MRS (magnetic resonance sounding) FID (frequency identity) signal noise inhibition method
CN105115622A (en) * 2015-08-12 2015-12-02 合肥工业大学 Denoising algorithm of fiber Raman temperature sensing system based on independent component analysis
CN105549097A (en) * 2015-12-22 2016-05-04 吉林大学 Transient electromagnetic signal power frequency and harmonic interference elimination method and apparatus thereof
CN105954170A (en) * 2016-05-04 2016-09-21 中国石油大学(华东) Nuclear magnetic porosity calculating method considering background noise signal quantity
CN106361328A (en) * 2016-10-21 2017-02-01 电子科技大学 Method for electroencephalogram signal extraction under magnetic resonance environment
CN106361328B (en) * 2016-10-21 2019-07-19 电子科技大学 A kind of EEG signals extracting method under magnetic resonance environment
CN106931905A (en) * 2017-03-09 2017-07-07 北京理工大学 A kind of digital Moiré patterns phase extraction method based on nonlinear optimization
CN106931905B (en) * 2017-03-09 2019-02-05 北京理工大学 A kind of digital Moiré patterns phase extraction method based on nonlinear optimization
CN106955112A (en) * 2017-03-17 2017-07-18 泉州装备制造研究所 Brain wave Emotion recognition method based on Quantum wavelet neural networks model
CN107092925A (en) * 2017-03-30 2017-08-25 中国人民解放军国防科学技术大学 Cerebral function magnetic resonance imaging blind source separation method based on SIM algorithms in groups
CN107092925B (en) * 2017-03-30 2019-09-17 中国人民解放军国防科学技术大学 Cerebral function magnetic resonance imaging blind source separation method based on SIM algorithm in groups
CN107121705A (en) * 2017-04-28 2017-09-01 中南大学 A kind of ground penetrating radar echo signals Denoising Algorithm compared based on automatic anti-phase correction and kurtosis value
CN108318764A (en) * 2018-03-28 2018-07-24 国网上海市电力公司 A kind of earthing or grounding means shock response test jamproof system and method
CN108846407A (en) * 2018-04-20 2018-11-20 太原理工大学 The nuclear magnetic resonance image classification method of brain network is not known based on independent element high order
CN108846407B (en) * 2018-04-20 2022-02-08 太原理工大学 Magnetic resonance image classification method based on independent component high-order uncertain brain network
CN109143389A (en) * 2018-08-01 2019-01-04 吉林大学 A kind of three-dimensional industrial frequency interference source for nuclear-magnetism quantifies orienting device and measurement method
CN110058312A (en) * 2018-10-22 2019-07-26 南方科技大学 A kind of method, apparatus and terminal device inhibiting the interference of earth magnetism near field noise
CN110058312B (en) * 2018-10-22 2020-07-31 南方科技大学 Method and device for inhibiting geomagnetic near-field noise interference and terminal equipment
CN109828318A (en) * 2019-01-25 2019-05-31 吉林大学 A kind of magnetic resonance depth measurement signal noise filtering method based on variation mode decomposition
CN111368778A (en) * 2020-03-13 2020-07-03 哈尔滨理工大学 Weak signal noise stripping method based on intelligent optimization algorithm
CN111368778B (en) * 2020-03-13 2022-05-27 哈尔滨理工大学 Weak signal noise stripping method based on intelligent optimization algorithm
CN112036453A (en) * 2020-08-14 2020-12-04 哈尔滨工程大学 Blind source separation method based on quantum rhinoceros search mechanism
CN112036453B (en) * 2020-08-14 2022-04-29 哈尔滨工程大学 Blind source separation method based on quantum rhinoceros search mechanism
CN112370066A (en) * 2020-09-30 2021-02-19 北京工业大学 Brain-computer interface method of stroke rehabilitation system based on generation of countermeasure network
CN112180454A (en) * 2020-10-29 2021-01-05 吉林大学 Magnetic resonance underground water detection random noise suppression method based on LDMM
CN112180454B (en) * 2020-10-29 2023-03-14 吉林大学 Magnetic resonance underground water detection random noise suppression method based on LDMM

Also Published As

Publication number Publication date
CN104459809B (en) 2017-04-12

Similar Documents

Publication Publication Date Title
CN104459809A (en) Full-wave nuclear magnetic resonance signal denoising method based on independent component analysis
US10677957B2 (en) Method for random noise reduction from MRS oscillating signal using joint algorithms of EMD and TFPF
CN107045149B (en) A kind of all-wave NMR signal noise filtering method based on double singular value decompositions
CN109828318B (en) Magnetic resonance sounding signal noise filtering method based on variational modal decomposition
Wu et al. Removal of multisource noise in airborne electromagnetic data based on deep learning
CN104777442B (en) A kind of nuclear magnetic resonance depth measurement FID signal noise suppressing method
Larsen Model-based subtraction of spikes from surface nuclear magnetic resonance data
CN108345039B (en) A method of eliminating adjacent frequency harmonic wave interference in ground nuclear magnetic resonance data
CN106646637A (en) Method for removing peak noise in nuclear magnetism signal
Ji et al. Noise reduction of grounded electrical source airborne transient electromagnetic data using an exponential fitting-adaptive Kalman filter
CN104614778A (en) Nuclear magnetic resonance underground water detection signal noise eliminating method based on independent component analysis (ICA)
CN107957566A (en) Magnetic resonance depth measurement method for extracting signal based on frequency selection singular spectrum analysis
Larsen et al. Processing of surface-nuclear magnetic resonance data from sites with high noise levels
CN109633761A (en) Magnetic resonance signal industrial frequency noise method for reducing based on wavelet modulus maxima method
CN104122590B (en) A kind of gas-oil detecting method based on electromagnetic survey and system
Li et al. Magnetotelluric signal-noise separation method based on SVM–CEEMDWT
CN103995293B (en) Method for detecting magnetic resonance sounding signals
Tian et al. Noise suppression method for magnetic resonance sounding signals based on double singular value decomposition
Jiang et al. Hydrocarbon detection based on empirical mode decomposition, teager-kaiser energy, and the cepstrum
Tian et al. Noise cancellation of a multi-reference full-wave magnetic resonance sounding signal based on a modified sigmoid variable step size least mean square algorithm
CN109885906B (en) Magnetic resonance sounding signal sparse noise elimination method based on particle swarm optimization
CN109782363B (en) Magnetic resonance signal denoising method based on time domain modeling and frequency domain symmetry
CN105891889A (en) Gravity abnormal boundary enhancement method and device
CN105675983A (en) Weak harmonic wave signal detection and reconstruction methods in strong chaotic background
CN106226077B (en) A kind of detection method of the periodical transient signal based on time-varying singular value decomposition

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20170412

Termination date: 20171030

CF01 Termination of patent right due to non-payment of annual fee