CN112099100B - Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering - Google Patents
Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering Download PDFInfo
- Publication number
- CN112099100B CN112099100B CN202010860696.0A CN202010860696A CN112099100B CN 112099100 B CN112099100 B CN 112099100B CN 202010860696 A CN202010860696 A CN 202010860696A CN 112099100 B CN112099100 B CN 112099100B
- Authority
- CN
- China
- Prior art keywords
- signal
- envelope
- equation
- imf
- mrs
- 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.)
- Active
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 38
- 238000000605 extraction Methods 0.000 title claims abstract description 21
- 230000003044 adaptive effect Effects 0.000 title claims description 17
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 25
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 22
- 238000001228 spectrum Methods 0.000 claims abstract description 7
- 238000000034 method Methods 0.000 claims description 52
- 230000008569 process Effects 0.000 claims description 31
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 11
- 238000005481 NMR spectroscopy Methods 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 3
- 238000009792 diffusion process Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 15
- 238000010586 diagram Methods 0.000 description 14
- 230000003595 spectral effect Effects 0.000 description 6
- 230000001629 suppression Effects 0.000 description 5
- 125000004435 hydrogen atom Chemical group [H]* 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 239000003673 groundwater Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- DDNCQMVWWZOMLN-IRLDBZIGSA-N Vinpocetine Chemical compound C1=CC=C2C(CCN3CCC4)=C5[C@@H]3[C@]4(CC)C=C(C(=O)OCC)N5C2=C1 DDNCQMVWWZOMLN-IRLDBZIGSA-N 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 230000005358 geomagnetic field Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 229960000744 vinpocetine Drugs 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V5/00—Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity
- G01V5/20—Detecting prohibited goods, e.g. weapons, explosives, hazardous substances, contraband or smuggled objects
- G01V5/22—Active interrogation, i.e. by irradiating objects or goods using external radiation sources, e.g. using gamma rays or cosmic rays
Landscapes
- Physics & Mathematics (AREA)
- High Energy & Nuclear Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
The invention belongs to the field of magnetic resonance sounding signal processing, in particular to a magnetic resonance sounding signal envelope extraction method based on self-adaptive local iterative filtering, which directly converts an MRS signal into an envelope signal f (t) by utilizing Hilbert transform and spectrum shift; determining a decomposition order of the envelope signal and a mask coefficient selected when a filtering interval is calculated by adopting a traversal mode for the envelope signal, respectively adopting an ALIF algorithm to extract IMF components through an inner loop and extract an effective MRS signal envelope through an outer loop to finish an extraction process for a real part and an imaginary part of the envelope signal f (t) according to the decomposition order and the mask coefficient, and extracting the effective MRS signal envelope through a margin combination mode. The problem of modal aliasing existing in EMD can be effectively avoided.
Description
Technical Field
The invention belongs to the field of Magnetic Resonance Sounding (MRS) signal processing, and particularly relates to a Magnetic Resonance Sounding signal envelope extraction method based on adaptive local iterative filtering.
Background
The surface nuclear Magnetic Resonance (MRS) is an advanced method capable of quantitatively determining the storage state of groundwater in the world at present because it directly acts on hydrogen protons in groundwater. Compared with other geophysical exploration methods, the MRS method can provide hydrogeological information such as the content and position of underground water and the porosity of a water storage layer medium without drilling. Therefore, in recent years, the method is widely applied to the aspects of underground water detection, evaluation, monitoring and early warning of hidden danger of geological disasters and the like. The method excites hydrogen protons in underground water through an artificially generated alternating electromagnetic field, and the excited hydrogen protons generate energy level transition. When the artificial field is removed, the acquisition system receives MRS signals generated when the hydrogen protons release energy to return, and detection of underground water is achieved. The general expression for a full-wave nuclear magnetic resonance signal is:
the equation contains four key parameters of the MRS signal: e0、
fL、
Respectively, an initial amplitude of the signal related to the water content of the formation, a relaxation time related to the porosity of the aquifer (the relaxation time of free water ranges from 30ms to 1000ms), larmor frequencies with different values for different zones, and an initial phase related to the conductivity of the aquifer.
Although the MRS method has the above advantages, in practical application, because the method uses a natural geomagnetic field as a background excitation field, shielding measures cannot be taken, so random noise, power frequency harmonic interference and accidental spike noise in the environment all affect a nano-volt level weak MRS signal, and effective extraction of the MRS signal is difficult to achieve.
In view of the above problems, a great deal of research work has been conducted by experts and scholars at home and abroad. The classical process flow comprises the following: firstly, removing spike noise by adopting a nonlinear energy operator or a statistical means; secondly, removing power line harmonics by adopting a harmonic modeling method based on a least square method theory; then, overlapping the data acquired for many times to suppress random noise; and finally, extracting envelope and gate integral processing of the superposed data, and providing the processed data to inversion software to realize the explanation of the underground hydrogeological information. For corrupted MRS signals, even those completely swamped by noise, suppression or removal of the noise is inevitably required to extract the required information. At present, for filtering power frequency harmonic waves and random noise: legchenko and Valla, in Removal of power-line harmonics from magnetic resonance measures (Journal of Applied Geophysics, Vol. 53: 103-. In the article "Adaptive noise cancellation of multichannel magnetic resonance signals" (geophysics Journal International, vol.191: 88-100.) by Dalgaard et al, it is proposed to use an Adaptive noise cancellation algorithm to process a multichannel MRS probe signal, so as to cancel out the correlated noise in the main channel and the reference channel. Tianbaofeng et al, in a paper "magnetic resonance signal noise suppression method based on reference coil and variable step size adaptation" (in the Proc. geophysical sciences, vol. 55, No. 7: 2462-. Ghanati et al, in the article "Joint application of a static optimization process empirical mode decomposition to MRS noise cancellation" (Journal of Applied geomatics, vol 111, 5: 110 + 120), proposed the extraction of effective MRS signal attenuation trend based on EMD and CEEMD methods. Larsen et al, in the paper "Noise cancel harmonic signal combining model-based removal of powerline harmonics and multichannel Wiener filtering" (geographic Journal International 2014, volume 196, 2: 828. sup. 836 page.), realizes effective Noise suppression by combining power frequency harmonic modeling with multichannel Wiener filtering, and proposes a method adopting delay harmonic modeling for the special case of same frequency harmonic interference. Qi Wang et al, in An paper "An alternative adaptive approach to handling co-frequency harmonic in surface nuclear magnetic resonance data" (Geophanic Journal International, 2018, Vol. 215, 1962 and page 1973.) propose that power frequency harmonic modeling based on the least squares theory is used to remove power line harmonic noise interference. In addition, some scholars propose to perform noise suppression on MRS signals by statistical means such as Singular Spectral Analysis (SSA) and bisingular value decomposition (D-SVD). The data processing means adopts various strategies to remove and suppress noise, and finally adopts an envelope extraction algorithm to extract the MRS signal attenuation envelope. Although the above algorithms all obtain significant effects under certain conditions, due to the difference of geographic environmental conditions, and the irregularity, instability and specific spatial distribution of noise sources, the application of the algorithms is limited, and a more effective and good algorithm needs to be continuously searched for and applied to the aspects of MRS signal noise elimination, suppression and information extraction.
Disclosure of Invention
Aiming at the problems that the complexity of noise contained in MRS signals and the conventional EMD algorithm easily generates mode aliasing, the invention provides a magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering, which can effectively avoid the mode aliasing problem of EMD.
The present invention is achieved in such a way that,
a magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering comprises the following steps:
step 1: a group of MRS signals F (t) acquired by a nuclear magnetic resonance sounding water detector are directly converted into envelope signals f (t) by Hilbert transform and spectrum shifting;
step 2: determining the decomposition order of the envelope signal in the step 1 and the mask coefficient selected when calculating the filtering interval by adopting a traversal mode for the envelope signal, wherein the mask coefficient is determined as a certain numerical value firstly, and then the initial amplitude E of the MRS signal is used0And MRS signal mean decay time
Determining an optimal decomposition order; after the resolution order is determined, the initial amplitude E of MRS signal is determined by traversing the coefficients between 1.0-3.0 at intervals of 0.20And MRS signal mean decay time
To determine the optimal mask coefficients;
and step 3: and (3) extracting an IMF component from the real part and the imaginary part of the envelope signal f (t) by adopting an ALIF algorithm through an inner loop and extracting an effective MRS signal envelope by ending an extraction process through an outer loop respectively according to the decomposition order and the mask coefficient selected in the step (2) in a margin combination mode.
Further, the inner loop process in step 3 comprises: calculating the convolution of an envelope signal f (t) to be decomposed and a filter function omega (t) to obtain a sliding operator gamma (f) (t), then extracting the solved sliding operator gamma (f) (t) from the envelope signal f (t), calculating to obtain a fluctuation operator kappa (f) (t), if the fluctuation operator kappa (f) (t) meets an IMF condition, taking kappa (f) (t) as an extracted IMF component c (t), and if the fluctuation operator kappa (f) (t) does not meet the IMF condition, continuously and repeatedly screening to obtain the IMF component of the fixed frequency component.
Further, in the step 3, the external circulation process includes removing the IMF components c (t) obtained in the internal circulation process from the envelope signals f (t) to be decomposed to obtain a margin γ (t), terminating iteration if the margin γ (t) presents a monotonous trend characteristic, otherwise, taking the γ (t) as an original signal, entering the internal circulation, continuing to repeat the above steps until the condition is met, terminating the circulation, and finally obtaining a plurality of IMF components.
Further, the inner circulation process in step 3 specifically includes:
step 3 a: selecting a value of an expected decomposition order M for a real part and an imaginary part and a value of a mask coefficient lambda for calculating a filtering interval;
and step 3 b: calculating a filtering interval l according to equation (1)n
In the formula, N represents the number of sampling points of envelope signal f (x), k represents the number of extreme points of signal to be decomposed, λ represents a scalar parameter, and λ ∈ [1, 3 ]]，
Represents rounding down;
and step 3 c: solving a filter function according to a Fokker-Planck equation of the formula (2)
Wherein p (x) and q (x) are both smoothly derivable functions and are in [ a, b ], wherein a < 0 < b, satisfying: q (a) q (b) 0, and q (x) > 0 for any x e (a, b);
P(a)＜0＜P(b)；
alpha and beta are called steady state coefficients, and the value range is (0, 1);
in the equation
Will produce a diffusion effect and drive the equation solution h (x) to move from the center point of the interval (a, b) to the end points a and b; at the same time, the user can select the desired position,
bringing h (x) from both ends of a and b toward the center of the interval (a, b), and when the two are balanced, the following are provided:
at this time, the equation has a non-zero solution h (x) and satisfies the condition:
for any x epsilon (a, b), h (x) is more than or equal to 0;
for any one
H (x) is 0;
representing that all equation solutions fall on the interval [ a, b ], solving the Fokker-Planck equation solution h (x) as a solved filter function omega (t), wherein the filter function omega (t) solves different analytic values along with the change of the interval [ a, b ], so that the self-adaptive solution of the ALIF algorithm to the filter function is realized;
and step 3 d: solving the sliding operator Γ (f) (x) using equation (4)
In the formula, ω (t) represents a filter function, lnRepresenting a filtering interval;
step 3 e: computing the fluctuation operator κ (f) (x) using equation (5)
κ(f)(x)＝f(x)-Γ(f)(x) (5)
And step 3 f: by using maleEquation (6) judges whether the fluctuation operator κ (f) (x) satisfies the condition of the IMF component, and if so, obtains the first IMF component c1(x) If not, taking kappa (f) (x) as the signal to be decomposed, and continuing to repeat the processes from the step 3b to the step 3e until the condition of IMF component is satisfied, stopping screening,
wherein, IMFs (x) represents the current IMF to be judged, IMFc (x) represents the previous IMF component, and δ represents the defined threshold, which ranges from 0.001 to 0.2.
Further, the external circulation process in the step 3 comprises the following steps:
step 3 g: obtaining the margin by equation (7)
γ1(x)＝f(x)-c1(x) (7)
Will gamma1(x) Repeating the inner loop process as the original signal to obtain the 2 nd IMF component, and using the formula (8) as the judgment gamman(x) The criterion for presenting a monotonic trend characteristic, when the value of epsilon is greater than 20, indicates gamman(x) If the trend characteristic is monotonous, the circulation is terminated to obtain a plurality of IMF components, gamman(x) Is the nth margin, c1(x) The first IMF component, where f (x) is represented by equation (9);
where f is the current original signal, γkIn general, the reasonable value of epsilon is 20;
wherein, ci(x) Is the ith IMF component.
Compared with the prior art, the invention has the beneficial effects that:
the invention selects a self-adaptive local iterative filtering algorithm as an extraction method of signal envelopes, the method comprises the steps of solving by utilizing the obtained filtering function and filtering interval to obtain a sliding operator, then extracting the sliding operator from an original signal, calculating to obtain a fluctuation operator, checking whether the fluctuation operator meets the condition of IMF component, continuously updating the filtering function and the filtering interval to obtain the IMF component until all IMF components are extracted from the original signal and the remaining margin has obvious monotone trend characteristics, and finally obtaining a plurality of intrinsic mode function components to realize effective decomposition of the signal. In view of the fact that the adaptive local iterative filtering has a solid mathematical foundation and can solve the problem of mode aliasing existing in the traditional EMD, the invention provides the method for extracting the envelope of the MRS signal by adopting the ALIF algorithm directly, and by means of the characteristic of algorithm self-adaptation, the filtering function is selected in a self-adaptive mode, so that the filtering process has smoothness, and artificial oscillation cannot be generated in the IMF component generated subsequently. In addition, in order to capture non-stationary changes in signal frequency and amplitude, the length of the filter is adjusted accordingly. Therefore, according to local conditions in the whole signal processing process, all signals are not considered together, unnecessary influence factors are reduced, and the accuracy and the reliability of a processing result are greatly improved. Therefore, the effective application of the method has profound theoretical significance for improving the anti-interference capability of the magnetic resonance sounding technology, has new application value for directly extracting the envelope of the magnetic resonance sounding signal, and has instructive effect on the thought inspiration of the theoretical method research in the field of signal processing.
Drawings
Fig. 1 is a block flow diagram of magnetic resonance sounding signal envelope extraction based on adaptive local iterative filtering;
FIG. 2 is a block flow diagram of an adaptive local iterative filtering algorithm;
fig. 3 shows an ideal MRS signal and frequency spectrum, fig. 3A is a time-amplitude diagram, and fig. 3B is a frequency-power spectral density diagram;
fig. 4 illustrates a noisy MRS signal and its frequency spectrum, fig. 4A is a time-amplitude diagram, and fig. 4B is a frequency-power spectral density diagram;
fig. 5 illustrates a full-wave signal of noisy MRS full-wave signal converted into an envelope signal, fig. 5A illustrates a time-amplitude diagram of a full-wave signal and its envelope signal, and fig. 5B illustrates a time-amplitude diagram of a real part and an imaginary part of an envelope signal.
Fig. 6 is an IMF component graph obtained by decomposition processing of the ALIF algorithm, fig. 6A is an IMF component time-amplitude graph obtained by decomposition of the real part of the original envelope signal, and fig. 6B is an IMF component time-amplitude graph obtained by decomposition of the imaginary part thereof.
Fig. 7 is a graph of the residual time-amplitude remaining after extraction of the IMF components, fig. 7A, 7B being real and imaginary, respectively.
Fig. 8 shows the signal envelope resulting from the residual combination of the real and imaginary parts, fig. 8A is a time-amplitude diagram, and fig. 8B is a frequency-amplitude diagram.
Fig. 9 shows the measured MRS signal and its power spectrum, fig. 9A is a time-amplitude graph, and fig. 9B is a frequency-power spectral density graph.
Fig. 10 measures an envelope signal finally obtained by the MRS signal.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is further described in detail with reference to the following embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
As shown in fig. 1, the magnetic resonance sounding signal envelope extraction algorithm based on adaptive local iterative filtering includes the following steps:
step 1: a group of MRS signals F (t) acquired by a nuclear magnetic resonance sounding water detector are directly converted into envelope signals f (t) by Hilbert transform and spectrum shifting;
step 2: determining the order of decomposition and the mask coefficient selected during the calculation of the filtering interval in a traversal manner, i.e. determining the mask coefficient as a certain value, and then determining the initial amplitude E of the MRS signal0And MRS signal mean decay time
Determining an optimal decomposition order; after the resolution order is determined, the initial amplitude E of MRS signal is determined by traversing the coefficients between 1.0-3.0 at intervals of 0.20And MRS signal mean decay time
To determine the optimal mask coefficients;
and step 3: and (3) according to the key parameters selected in the step (2), including the decomposed order and the mask coefficient selected when the filtering interval is calculated, respectively processing the real part and the imaginary part by adopting an ALIF algorithm, and extracting a good and effective MRS signal envelope in a margin combination mode.
wherein the internal circulation process is as follows: solving to obtain a sliding operator gamma (f) (t) by calculating the convolution of the envelope signal f (t) to be decomposed and the filter function omega (t), then extracting the solved sliding operator gamma (f) (t) from the envelope signal f (t), calculating to obtain a fluctuation operator kappa (f) (t), and if the fluctuation operator kappa (f) (t) meets the IMF condition, obtaining the extracted IMF component kappa (f) (t). Under the general condition, the fluctuation operators obtained by the first calculation do not meet the condition of the IMF component, and the IMF component of the fixed frequency component can be obtained by continuously and repeatedly screening in a circulating way. Therefore, the inner loop is stopped until the condition of the IMF component is satisfied.
The external circulation process is as follows: and (3) removing the IMF components c (t) obtained in the inner loop process from the signals f (t) to be decomposed to obtain the residual gamma (t), terminating iteration if the gamma (t) presents an obvious monotonous trend characteristic, otherwise, taking the gamma (t) as an original signal, entering the inner loop, continuously repeating the steps until the conditions are met, terminating the loop, and finally obtaining a plurality of IMF components.
As shown in fig. 2, the Adaptive Local Iterative Filtering (ALIF) algorithm provided by the present invention performs processing on a noisy MRS signal, that is, the specific steps of the inner loop process and the outer loop process are as follows:
step 3 a: selecting a value of an expected decomposition order M for a real part and an imaginary part and a value of a mask coefficient lambda for calculating a filtering interval;
and step 3 b: calculating a filtering interval l according to equation (1)n
In the formula, N represents the number of sampling points of envelope signal f (x), k represents the number of extreme points of signal to be decomposed, λ represents a scalar parameter, and λ ∈ [1, 3 ]]，
Indicating a rounding down.
And step 3 c: solving a filter function according to a Fokker-Planck equation of the formula (2)
Wherein p (x) and q (x) are both smoothly derivable functions and satisfy, at [ a, b ] (where a < 0 < b):
(1) q (a) q (b) 0, and q (x) > 0 for any x e (a, b);
(2)p(a)＜0＜p(b)。
alpha and beta are called steady state coefficients, and the value range is (0, 1).
In the equation
Will produce a diffusion effect and drive the equation solution h (x) to move from the center point of the interval (a, b) to the end points a and b; at the same time, the user can select the desired position,
h (x) from both ends of a and b to the center of the interval (a, b). When the two are in equilibrium, there are:
at this time, the equation has a non-zero solution h (x) and satisfies the condition:
(1) for any x epsilon (a, b), h (x) is more than or equal to 0;
(2) for any one
H (x) is 0;
this means that the equation solution set falls on the interval [ a, b ] completely, and the solution h (x) of the Fokker-Planck equation is the sought filter function ω (t). And omega (t) is solved to obtain different analytic values along with the change of the intervals (a, b), so that the adaptive solution of the ALIF algorithm to the filter function is realized.
And step 3 d: solving the sliding operator Γ (f) (x) using equation (4)
In the formula, ω (t) represents a filter function, lnRepresenting a filtering interval;
step 3 e: computing the fluctuation operator κ (f) (x) using equation (5)
κ(f)(x)＝f(x)-Γ(f)(x) (5)
And step 3 f: judging whether the fluctuation operator kappa (f) (x) meets the condition of the IMF component by using the formula (6), and if so, obtaining a first IMF component c1(x) If not, taking kappa (f) (x) as a signal to be decomposed, and continuing to repeat the processes from the step 3b to the step 3e until the condition of IMF component is met, stopping screening.
In the formula, IMFs (x) represents the current IMF to be determined, and IMFc (x) represents the previous IMF component. δ represents a defined threshold, ranging between 0.001 and 0.2.
Step 3 g: step 3 g: obtaining the margin by equation (7)
γ1(x)＝f(x)-c1(x) (7)
Will gamma1(x) Repeating the inner loop process as the original signal to obtain the 2 nd IMF component, and using the formula (8) as the judgment gamman(x) The criterion for presenting a monotonic trend characteristic, when the value of epsilon is greater than 20, indicates gamman(x) If the trend characteristic is monotonous, the circulation is terminated to obtain a plurality of IMF components, gamman(x) Is the nth margin, c1(x) The first IMF component, where f (x) is represented by equation (9);
where f is the current original signal, γkIn general, the reasonable value of epsilon is 20;
wherein, ci(x) Is the ith IMF component.
Example 1
This example is a simulation experiment of the method of the present invention conducted in the MATLAB 7.0 programming environment.
A magnetic resonance sounding signal envelope extraction algorithm simulation experiment based on adaptive local iterative filtering comprises the following steps:
step (1): by using
Constructing Larmor frequency of 2114Hz and amplitude E0180nV, relaxation time
The ideal MRS signal is 100ms, as shown in fig. 3, where fig. 3A is a time-amplitude graph and fig. 3B is a frequency-power spectral density graph. In addition, Gaussian random white noise with the amplitude of 10nV is added, and the amplitude isA noisy MRS signal generated by 80 power frequency harmonics with phases (-pi, pi) randomly selected within (0, 40) nV and fundamental frequencies between (49.9, 50.1) is shown in fig. 4, where fig. 4A is a time-amplitude graph and fig. 4B is a frequency-power spectral density graph;
step (2): first, the generated noise-containing MRS full-wave signal is converted into an envelope signal, as shown in fig. 5, fig. 5A is a time-amplitude diagram of the full-wave signal and its envelope signal, and fig. 5B is a time-amplitude diagram of the real part and imaginary part of the envelope signal. Then, the decomposition process of the ALIF algorithm is performed to obtain 7 IMF components, as shown in fig. 6, and fig. 6A, B shows an exploded view of the real part and the imaginary part of the full-wave signal after being converted into the envelope signal, respectively. Finally, after extracting the IMF component from the original signal, the remaining margin left is the target to be acquired, as shown in fig. 7, and fig. 7A and 7B are the margin of the real part and the margin of the imaginary part, respectively.
And (3): through the above process, further, the signal envelope can be directly obtained by combining the residuals of the real part and the imaginary part, as shown in fig. 8, where fig. 8A is a time-amplitude diagram and fig. 8B is a frequency-amplitude diagram.
In order to verify the practicability of the method, parameter extraction is carried out on the obtained MRS signal envelope, and the initial amplitude and the relaxation time of the obtained key parameters are respectively E0＝173.2983nV，T2 *The relative error is 3.7231% and 9.3126% respectively, which satisfy the application requirement 109.3 ms.
Example 2
In this embodiment, the MRS signal collected in the field of the cultural square in the vinpocetine city is used as the processing target of the method of the present invention.
The magnetic resonance sounding signal envelope extraction algorithm based on the adaptive local iterative filtering, referring to fig. 1, comprises the following steps:
step (1): a set of original noisy MRS signals acquired by a magnetic resonance depth sounding (MRS) water detector, the local larmor frequency is 2352Hz, as shown in fig. 9, fig. 9A is a time-amplitude graph, and fig. 9B is a frequency-power spectral density graph.
Step (2): and converting the acquired original full-wave MRS signal into an envelope signal, then carrying out decomposition processing of an ALIF algorithm, wherein the decomposition order is 7, and finally obtaining the remaining margin.
And (3): after the above processing procedure, the MRS signal envelope may be directly extracted, as shown in fig. 10.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.
Claims (2)
1. A magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering is characterized by comprising the following steps:
step 1: a group of MRS signals F (t) acquired by a nuclear magnetic resonance sounding water detector are directly converted into envelope signals f (t) by Hilbert transform and spectrum shifting;
step 2: determining the decomposition order of the envelope signal in the step 1 and the mask coefficient selected when calculating the filtering interval by adopting a traversal mode for the envelope signal, wherein the mask coefficient is determined as a certain numerical value firstly, and then the initial amplitude E of the MRS signal is used0And MRS signal mean decay time
Determining an optimal decomposition order; after the resolution order is determined, the initial amplitude E of MRS signal is determined by traversing the coefficients between 1.0-3.0 at intervals of 0.20And MRS signal mean decay time
To determine the optimal mask coefficients;
and step 3: according to the decomposition order and the mask coefficient selected in the step 2, respectively adopting an ALIF algorithm to extract IMF components through inner circulation and extract the process through outer circulation termination, and extracting effective MRS signal envelope through a margin combination mode for the real part and the imaginary part of the envelope signal f (t);
the internal circulation process in the step 3 comprises the following steps: solving to obtain a sliding operator gamma (f) (t) by calculating the convolution of an envelope signal f (t) to be decomposed and a filter function omega (t), then extracting the solved sliding operator gamma (f) (t) from the envelope signal f (t), calculating to obtain a fluctuation operator kappa (f) (t), if the fluctuation operator kappa (f) (t) meets an IMF condition, taking kappa (f) (t) as an extracted IMF component c (t), and if the fluctuation operator kappa (f) (t) does not meet the IMF condition, continuously and repeatedly screening to obtain an IMF component of a fixed frequency component;
in the step 3, the outer loop process comprises the steps of removing IMF components c (t) obtained in the inner loop process from envelope signals f (t) to be decomposed to obtain a margin gamma (t), terminating iteration if the margin gamma (t) presents a monotonous trend characteristic, otherwise, taking the gamma (t) as an original signal, entering the inner loop, continuously repeating the steps until the conditions are met, terminating the loop, and finally obtaining a plurality of IMF components;
the internal circulation process in the step 3 specifically comprises the following steps:
step 3 a: selecting a value of an expected decomposition order M for a real part and an imaginary part and a value of a mask coefficient lambda for calculating a filtering interval;
and step 3 b: calculating a filtering interval l according to equation (1)n
In the formula, N represents the number of sampling points of envelope signal f (x), k represents the number of extreme points of signal to be decomposed, λ represents a scalar parameter, and λ ∈ [1, 3 ]]，
Represents rounding down;
and step 3 c: solving a filter function according to a Fokker-Planck equation of the formula (2)
Wherein p (x) and q (x) are both smoothly derivable functions and are in [ a, b ], wherein a < 0 < b, satisfying: q (a) q (b) 0, and q (x) > 0 for any x e (a, b);
p(a)＜0＜p(b)；
alpha and beta are called steady state coefficients, and the value range is (0, 1);
in the equation
Will produce a diffusion effect and drive the equation solution h (x) to move from the center point of the interval (a, b) to the end points a and b; at the same time, the user can select the desired position,
bringing h (x) from both ends of a and b toward the center of the interval (a, b), and when the two are balanced, the following are provided:
at this time, the equation has a non-zero solution h (x) and satisfies the condition:
for any x epsilon (a, b), h (x) is more than or equal to 0;
for any one
H (x) is 0;
representing that all equation solutions fall on the interval [ a, b ], solving the Fokker-Planck equation solution h (x) as a solved filter function omega (t), wherein the filter function omega (t) solves different analytic values along with the change of the interval [ a, b ], so that the self-adaptive solution of the ALIF algorithm to the filter function is realized;
and step 3 d: solving the sliding operator Γ (f) (x) using equation (4)
In the formula, ω (t) represents a filter function, lnRepresenting a filtering interval;
step 3 e: computing the fluctuation operator κ (f) (x) using equation (5)
κ(f)(x)＝f(x)-Γ(f)(x) (5)
And step 3 f: judging whether the fluctuation operator kappa (f) (x) meets the condition of the IMF component by using the formula (6), and if so, obtaining a first IMF component c1(x) If not, taking kappa (f) (x) as the signal to be decomposed, and continuing to repeat the processes from the step 3b to the step 3e until the condition of IMF component is satisfied, stopping screening,
wherein, IMFs (x) represents the current IMF to be judged, IMFc (x) represents the previous IMF component, and δ represents the defined threshold, which ranges from 0.001 to 0.2.
2. The method of claim 1, wherein the outer loop process in step 3 comprises:
step 3 g: obtaining the margin by equation (7)
γ1(x)＝f(x)-c1(x) (7)
Will gamma1(x) Repeating the inner loop process as the original signal to obtain the 2 nd IMF component, and using the formula (8) as the judgment gamman(x) The criterion for presenting a monotonic trend characteristic, when the value of epsilon is greater than 20, indicates gamman(x) If the trend characteristic is monotonous, the circulation is terminated to obtain a plurality of IMF components, gamman(x) Is the nth margin, c1(x) The first IMF component, where f (x) is represented by equation (9);
where f is the current original signal, γkThe value of epsilon is 20 for the currently obtained margin;
wherein, ci(x) Is the ith IMF component.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010860696.0A CN112099100B (en) | 2020-08-25 | 2020-08-25 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010860696.0A CN112099100B (en) | 2020-08-25 | 2020-08-25 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112099100A CN112099100A (en) | 2020-12-18 |
| CN112099100B true CN112099100B (en) | 2021-08-27 |
Family
ID=73754468
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010860696.0A Active CN112099100B (en) | 2020-08-25 | 2020-08-25 | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112099100B (en) |
Family Cites Families (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102368095B (en) * | 2011-09-10 | 2013-03-27 | 吉林大学 | Extraction method for relaxation time spectrum of nuclear magnetic resonance detection signal for underground water by utilizing multi exponent fitting technology |
| FR3016966A1 (en) * | 2014-01-29 | 2015-07-31 | Centre Nat Rech Scient | DEVICE AND METHOD FOR ULTRASONIC CONTROL OF BONE QUALITY |
| CN106198018B (en) * | 2016-06-29 | 2018-05-25 | 潍坊学院 | A kind of EEMD of rotating machinery and smooth iteration envelope Analysis Method |
| CN106198013B (en) * | 2016-06-29 | 2018-05-29 | 潍坊学院 | A kind of envelope Analysis Method based on empirical mode decomposition filtering |
| CN107783200B (en) * | 2017-11-21 | 2019-06-07 | 吉林大学 | A kind of all-wave magnetic resonance signal random noise method for reducing for combining EMD and TFPF algorithm |
| CN108459353B (en) * | 2018-03-26 | 2019-06-21 | 吉林大学 | Faint magnetic resonance signal extracting method and device under a kind of electromagnetic noise background |
| CN109470135A (en) * | 2018-11-12 | 2019-03-15 | 吉林大学 | CSAMT data inactivity bearing calibration |
| CN109827697A (en) * | 2019-03-19 | 2019-05-31 | 东南大学 | Suspension cable time-varying Suo Li recognition methods based on local mean value mode decomposition |
| CN111239837B (en) * | 2020-02-20 | 2021-01-05 | 吉林大学 | Ground magnetic resonance signal parameter extraction method based on MCMC |
| CN111352163B (en) * | 2020-03-03 | 2021-04-02 | 吉林大学 | Magnetotelluric depth sounding-based static effect correction method and system |
-
2020
- 2020-08-25 CN CN202010860696.0A patent/CN112099100B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN112099100A (en) | 2020-12-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102998706B (en) | Method and system for attenuating seismic data random noise | |
| CN109828318B (en) | Magnetic resonance sounding signal noise filtering method based on variational modal decomposition | |
| Liu et al. | Seismic data interpolation beyond aliasing using regularized nonstationary autoregression | |
| CN104777442B (en) | A kind of nuclear magnetic resonance depth measurement FID signal noise suppressing method | |
| CN107783200A (en) | Joint EMD and TFPF algorithms a kind of all-wave magnetic resonance signal random noise method for reducing | |
| CN106772646B (en) | A kind of ground nuclear magnetic resonance method for extracting signal | |
| CN110388570B (en) | VMD-based self-adaptive noise reduction method and application thereof in water supply pipeline leakage positioning | |
| CN110095773A (en) | The two-parameter inversion method of the multiple dimensioned Full wave shape of Ground Penetrating Radar | |
| CN109633761B (en) | Magnetic resonance signal power frequency noise reduction method based on wavelet transformation modulus maximum value method | |
| CN108345039B (en) | A method of eliminating adjacent frequency harmonic wave interference in ground nuclear magnetic resonance data | |
| CN109885903B (en) | Model-based ground nuclear magnetic resonance signal spike noise removing method | |
| CN112882115B (en) | Magnetotelluric signal denoising method and system based on GWO optimized wavelet threshold | |
| CN113568058B (en) | Magnetotelluric signal-noise separation method and system based on multi-resolution singular value decomposition | |
| CN114091538B (en) | Intelligent noise reduction method for discrimination loss convolutional neural network based on signal characteristics | |
| CN108828658A (en) | A kind of ocean bottom seismic data reconstructing method | |
| Shi et al. | Signal extraction using complementary ensemble empirical mode in pipeline magnetic flux leakage nondestructive evaluation | |
| CN114970646A (en) | Artificial source electromagnetic pseudorandom signal detrending and noise identification method | |
| CN112099100B (en) | Magnetic resonance sounding signal envelope extraction method based on adaptive local iterative filtering | |
| CN110703089B (en) | Wavelet threshold denoising method for low-frequency oscillation Prony analysis | |
| Huang et al. | Frequency–space-dependent smoothing regularized nonstationary predictive filtering | |
| CN115826068A (en) | MRS signal envelope extraction method based on self-adaptive Gaussian filtering | |
| CN116401513A (en) | Magnetic resonance power frequency harmonic noise suppression method based on depth residual error network | |
| Lin et al. | Removal of a series of spikes from magnetic resonance sounding signal by combining empirical mode decomposition and wavelet thresholding | |
| Ma et al. | A Global and Multi-Scale Denoising Method Based on Generative Adversarial Network for DAS VSP Data | |
| Xia et al. | Automatic demarcation of sequence stratigraphy using the method of well logging multiscale data fusion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |