CN107957566B - Magnetic resonance depth measurement method for extracting signal based on frequency selection singular spectrum analysis - Google Patents

Abstract

The present invention is a kind of magnetic resonance depth measurement method for extracting signal of singular spectrum analysis based on frequency selection.MRS signal regional known to water instrument acquisition Larmor frequency, which is visited, first with nuclear magnetic resonance depth measurement finds the position of the corresponding MRS signal of Larmor frequency on power spectrum based on power spectrumanalysis after carrying out partial noise inhibition by broadband bandpass filter；Then, the singular spectrum analysis selected based on frequency is carried out to extract MRS signal.Singular spectrum analysis based on frequency selection includes being embedded in, RSVD decomposition, carrying out matrix reconstruction according to the corresponding singular value of MRS signal amplitude selection and diagonally equalize four steps.The present invention is able to solve effectively filtering out for random noise in noisy MRS signal, spike noise and industrial frequency harmonic interference, realize that complicated very noisy interferes effective extraction of lower MRS signal, compared with traditional MRS signal antinoise method, have many advantages, such as that arithmetic speed is fast, signal-to-noise ratio is high, practical.

Description

Magnetic Resonance Sounding Signal Extraction Method Based on Frequency Selective Singular Spectrum Analysis

technical field

The invention relates to the technical field of magnetic resonance sounding (also known as ground nuclear magnetic resonance) groundwater detection signal noise filtering and parameter extraction, in particular to a method for extracting magnetic resonance sounding signals based on frequency-selective singular spectrum analysis.

Background technique

Magnetic Resonance Sounding (MRS) technology is the only method to directly and quantitatively characterize water content and pore structure. Its basic principle is to detect groundwater by detecting the MRS signal response generated by the hydrogen-proton resonance transition in groundwater. The order of magnitude of the MRS signal is nanovolts, which is very weak, and is easily interfered by random noise, power frequency harmonics and spikes in the environment, which will affect the quality of the detected MRS signal, which in turn will affect the extraction of the characteristic parameters of the MRS signal. The accuracy of the inversion results is reduced, resulting in inaccurate assessment of the regional water resource content and storage medium composition.

At present, experts and scholars at home and abroad have carried out a lot of research work on the theory and method of noise filtering and signal extraction of magnetic resonance sounding signals. In "Symmetry based frequency domain processing to remove harmonic noise from surface nuclear magnetic resonance measurements"("Geophysical Journal International", No. 208, 2017: 724-736), Annette Hein et al. proposed a method based on the symmetry of the MRS signal in the frequency domain. characteristics to eliminate the method of power frequency noise, but can not deal with spike noise and random noise; Journal of Central South University of Technology "2017 No. 24: 900-911) proposed a multi-channel full-wave NMR signal reference cancellation method based on the improved s-function variable step size LMS, but when the main channel is related to the reference channel When the performance is poor, the treatment effect is not good. In "Magnetic Resonance Sounding Signal Detection Based on Bayes Bootstrap Statistical Noise Reduction Method"("Journal of Central South University (Natural Science Edition)" 2014, Vol. 45, No. 9: 3144-3149), Zhang Hairu et al. Two optimal estimation points were extracted from the MRS signal components to reconstruct the MRS signal, but the error of the signal parameter amplitude e ₀ and the average relaxation time T ₂ ^* extracted after processing was relatively large.

CN106772646A discloses "a ground nuclear magnetic resonance signal extraction method", which can adaptively find power frequency harmonics with a fundamental frequency of 49.9Hz to 50.1Hz, obtain the autocorrelation expression of the MRS signal, and quickly and effectively realize Separation of signal and noise, but this method is only aimed at dealing with power frequency harmonic noise in nuclear magnetic resonance signals; CN104459809A discloses "a method for filtering noise of full-wave nuclear magnetic resonance signals based on independent component analysis", using independent component analysis algorithm The power frequency noise is eliminated, and the digital quadrature method is used to construct a virtual input channel signal to solve the problem of underdetermined blind source separation, but this method is powerless against strong random noise and spike noise interference. CN106226407A discloses "an online preprocessing method for ultrasonic echo signals based on singular spectrum analysis", which uses singular spectrum analysis for echo signal preprocessing in ultrasonic online detection technology, and can automatically realize echo signal processing in ultrasonic detection. Noise removal and separation and extraction of signal components in different frequency bands; CN106404386A discloses "a method for collecting, extracting and diagnosing early fault characteristic signals of gearboxes", which uses singular spectrum analysis for fault diagnosis. It can be seen that singular spectrum analysis has been successfully applied to various fields of signal processing, but it has not been applied to noise filtering of MRS signals.

Contents of the invention

The technical problem to be solved by the present invention is to provide a method for extracting magnetic resonance sounding signals based on frequency-selective singular spectrum analysis to solve the problem of MRS signal extraction caused by the influence of complex peak noise, power frequency harmonics and random noise in the environment. .

The present invention is achieved in this way, based on the magnetic resonance sounding signal extraction method of frequency selected singular spectrum analysis (Frequency chosen of Singular-Spectrum-Analysis, F-SSA), the method comprises the following steps:

Step (1): Collect a set of observed MRS signals X ₁ (t) with a known Larmor frequency using a nuclear magnetic resonance sounding (MRS) water detector;

Step (2): pass the collected observed MRS signal X ₁ (t) through a broadband band-pass filter to obtain X ₂ (t);

Step (3): Perform power spectrum analysis on the X ₂ (t) obtained through the broadband bandpass filter, arrange the amplitudes of each frequency in descending order, and find the sorting position of the signal amplitude corresponding to the Larmor frequency;

Step (4): Perform singular spectrum analysis on X ₂ (t) obtained through a broadband bandpass filter to extract the MRS signal. The singular spectrum analysis includes: embedding to obtain the trajectory matrix H; performing RSVD decomposition on the trajectory matrix H to obtain Singular values and singular vector U arranged in descending order; according to the sorting position obtained in step (3), select two corresponding singular values for matrix reconstruction to obtain matrix C; diagonally average matrix C to obtain extracted MRS signal Y(t);

Further, the specific steps of embedding described in step (4) are:

The signal X ₂ (t) is a one-dimensional real sequence of length N, X ₂ (t)=(x ₁ ,x ₂ ,...,x _N ), the positive integer L is the length of the sliding window, 1<L<N , the value of L is determined by the following formula

By embedding the original sequence signal X ₂ (t) to form P vectors, each vector can be represented by h _i

h _i ＝(x _i ,x _i+1 ,…,x _i+L-1 ) ^T

Where P=N-L+1, i=1,2,...,P, the result of the mapping forms a trajectory matrix H:

Further, the concrete steps of the RSVD decomposition described in step (4) are:

1), set parameter k and parameter w, k is the k singular values taken for approximating the reconstruction matrix, w is a parameter used to ensure that the conditions of the reconstruction matrix are established, where w>k, w<L, w <P;

2), construct a Gaussian random matrix G _L×w with 0 mean and 1 variance;

3), Calculate the sampling matrix M _P×w of the trajectory matrix H _L×P

4) Select the first k singular values of the sampling matrix M _P×w and perform SVD decomposition to obtain the left singular vector Q _P×k of the orthogonal matrix, the diagonal matrix Z _k×k , and the right singular vector

5), construct matrix T _L×k :

T _L×k ＝H _L×P Q _P×k

6) Perform SVD decomposition on matrix T _L×k to obtain its left singular vector U _L×L , diagonal matrix Σ _L×k and right singular vector

7) Calculation matrix V _P×k :

V _{P × k} = Q _{P × k} O _{k × k} ;

8), calculate the approximate reduced rank matrix of H _L×P

Use RSVD decomposition to obtain the approximate reduced rank matrix of H, and the obtained approximate reduced rank matrix must meet the following conditions:

λ _k+1 is the k+1th singular value of matrix H;

The diagonal matrix Σ _L×k is a matrix of order L×k, and its main diagonal elements are the k singular values of the approximate reduced-rank matrix of the MRS signal; the k singular values are arranged in descending order (λ ₁ ≥λ ₂ ≥... ≥λ _k ), the corresponding orthogonal singular vector set U′=(u ₁ ,u ₂ ,...u _k ), the contribution rate of each singular value is:

Draw a singular spectrum according to the above formula.

Further, according to the sorting positions described in step (4), the specific steps for selecting two corresponding singular values to carry out matrix reconstruction are:

According to step (3), find the sorting v of the magnitude of the MRS signal amplitude, and select the 2v-1 and 2v singular values in the singular spectrum;

According to the selected singular value reconstruction matrix, the reconstruction steps are as follows:

1), calculate the right singular vector W of the reconstructed signal matrix:

Among them, j=2v-1,2v;

2), reconstruct signal matrix C

Further, the specific steps of diagonal averaging in the described step (4) are:

Singular spectrum analysis reconstructed signal matrix {c ₁ ,c ₂ ,…,c ^q ,…,c ^P } is transformed into corresponding reconstructed sequence {g ¹ ,g ² ,…,g ^q ,… by diagonal averaging ,g ^P }, where the sequence g ^q represents the qth singular spectrum analysis reconstruction sequence, The process is as follows:

Among them, c ^q represents the qth singular spectrum analysis reconstructed signal matrix, is the element in row m and column n of matrix c ^q , Represents the dth element of the reconstructed sequence g ^q , q=1,2,...,P, d=1,2,...,N; the reconstructed sequence {g ₁ ,g ₂ ,...,g ^P } accumulate and sum to obtain the denoised MRS signal sequence Where Y(t)={y(t ₁ ), y(t ₂ ), . . . , y(t _N )}={y ₁ , y ₂ , . . . , y _N }.

Compared with the prior art, the present invention has the beneficial effect that the present invention proposes a method for extracting magnetic resonance sounding signals based on frequency-selective singular spectrum analysis. For the detection data collected by a single channel, by finding the MRS signal amplitude in the power spectrum Sorting the position in, using the corresponding singular value for matrix reconstruction and signal recovery, can remove the interference of spike noise, power frequency harmonic and random noise at one time, and realize the effective extraction of MRS signal. The method of the invention solves the difficult problem that the signals caused by peak noise, power frequency harmonics and random noise are difficult to be effectively extracted in the work of magnetic resonance sounding to find water. Limitation saves a lot of financial and material resources, and opens up a new world of singular spectrum analysis in the field of NMR signal denoising.

Description of drawings

Fig. 1 is the process block diagram of the extraction method of the magnetic resonance sounding signal based on the singular spectrum analysis of frequency selection of the present invention;

Fig. 2 is a flowchart block diagram of the RSVD decomposition algorithm of the present invention;

Fig. 3 is time domain and power spectrum thereof of noisy MRS signal and pure MRS signal of the present invention, wherein (a) is time domain spectrum, (b) is power spectrum;

Fig. 4 is the time domain and its power spectrum before and after processing by the noise-containing MRS signal bandpass filter of the present invention, wherein (a) is the time domain spectrum, and (b) is the power spectrum;

Fig. 5 selects band-pass filter curve for the example 1 of the present invention;

Fig. 6 is the singular spectrogram of example 1 of the present invention;

Fig. 7 is the time domain and its power spectrum before and after processing the simulated MRS signal F-SSA of the present invention, wherein (a) is the time domain spectrum, (b) is the power spectrum;

Fig. 8 is the time domain and its power spectrum before and after the MRS signal band-pass filter processing of the actual measurement of the present invention, wherein (a) is the time domain spectrum, (b) is the power spectrum;

Fig. 9 selects the band-pass filter curve for the example of the present invention 2;

Fig. 10 is the singular spectrogram of Example 2 of the present invention;

Fig. 11 is the time domain and power spectrum of the measured MRS signal before and after F-SSA processing in the present invention, wherein (a) is the time domain spectrum, and (b) is the power spectrum.

Detailed ways

In order to make the object, technical solution and advantages of the present invention more clear, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention.

As shown in Figure 1, the method for extracting magnetic resonance sounding signals based on frequency-selective singular spectrum analysis includes the following steps:

Step (1): Collect a set of observed MRS signals X ₁ (t) with a known Larmor frequency using a nuclear magnetic resonance sounding (MRS) water detector;

Step (2): pass the collected observed MRS signal X ₁ (t) through a broadband band-pass filter to obtain X ₂ (t);

Step (3): Perform power spectrum analysis on the X ₂ (t) obtained through the broadband bandpass filter, arrange the amplitudes of each frequency in descending order, and find the sorting position of the signal amplitude corresponding to the Larmor frequency;

Step (4): performing singular spectrum analysis on the X ₂ (t) obtained through the broadband band-pass filter to extract the MRS signal. The singular spectrum analysis includes: embedding to obtain the trajectory matrix H; decomposing the trajectory matrix H by RSVD to obtain the singular values and singular vector U in descending order; according to the sorting position obtained in step (3), select two corresponding singular Values are matrix reconstructed to obtain matrix C; matrix C is diagonally averaged to obtain the extracted MRS signal Y(t);

The specific steps of embedding described in a magnetic resonance sounding signal extraction method based on frequency-selective singular spectrum analysis are:

The signal X ₂ (t) is a one-dimensional real sequence of length N, X ₂ (t)=(x ₁ ,x ₂ ,...,x _N ), the positive integer L is the length of the sliding window, 1<L<N , the value of L is determined by the following formula

By embedding the original sequence signal X ₂ (t) to form P vectors, each vector can be represented by h _i

h _i ＝(x _i ,x _i+1 ,…,x _i+L-1 ) ^T

Where P=N-L+1, i=1,2,...,P, the result of the mapping forms a trajectory matrix H:

As shown in Figure 2, the specific steps of the RSVD decomposition provided by the present invention are:

1), set parameter k and parameter w, k is the k singular values taken for approximating the reconstruction matrix, w is a parameter used to ensure that the conditions of the reconstruction matrix are established, where w>k, w<L, w <P;

2), construct a Gaussian random matrix G _L×w with 0 mean and 1 variance;

3), Calculate the sampling matrix M _P×w of the trajectory matrix H _L×P

4) Select the first k singular values of the sampling matrix M _P×w and perform SVD decomposition to obtain the left singular vector Q _P×k of the orthogonal matrix, the diagonal matrix Z _k×k , and the right singular vector

5), construct matrix T _L×k :

T _L×k ＝H _L×P Q _P×k

6) Perform SVD decomposition on matrix T _L×k to obtain its left singular vector U _L×L , diagonal matrix Σ _L×k and right singular vector

7) Calculation matrix V _P×k :

V _{P × k} = Q _{P × k} O _{k × k} ;

8) Calculate the approximate rank-reduced matrix H～L× _P of H _L×P :

Use RSVD decomposition to obtain the approximate reduced rank matrix of H, and the obtained approximate reduced rank matrix must meet the following conditions:

λ _k+1 is the (k+1)th singular value of matrix H;

The diagonal matrix Σ _L×k is a matrix of order L×k, and its main diagonal elements are the k singular values of the approximate reduced-rank matrix of the MRS signal; the k singular values are arranged in descending order (λ ₁ ≥λ ₂ ≥... ≥λ _k ), the corresponding orthogonal singular vector set U′=(u ₁ ,u ₂ ,...u _k ), the contribution rate of each singular value is:

Draw a singular spectrum according to the above formula.

The specific steps of matrix reconstruction described in a method for extracting full-wave nuclear magnetic resonance signals based on frequency-selective singular spectrum analysis are:

According to step (3), find the sorting v of the magnitude of the MRS signal amplitude, and select the 2v-1 and 2v singular values in the singular spectrum;

According to the selected singular value reconstruction matrix, the reconstruction steps are as follows:

1) Calculate the right singular vector W of the reconstructed signal matrix

Among them, j=2v-1,2v;

2), reconstruct signal matrix C

The specific steps of diagonal averaging described in a full-wave nuclear magnetic resonance signal extraction method based on frequency-selective singular spectrum analysis are:

Singular spectrum analysis reconstructed signal matrix {c ₁ ,c ₂ ,…,c ^q ,…,c ^P } is transformed into corresponding reconstructed sequence {g ¹ ,g ² ,…,g ^q ,… by diagonal averaging ,g ^P }, where the sequence g ^q represents the qth singular spectrum analysis reconstruction sequence, The process is as follows:

Among them, c ^q represents the qth singular spectrum analysis reconstructed signal matrix, is the element in row m and column n of matrix c ^q , Represents the dth element of the reconstructed sequence g ^q , q=1,2,...,P, d=1,2,...,N. Accumulate and sum the reconstructed sequence {g ₁ ,g ₂ ,…,g ^P } to obtain the denoised MRS signal sequence Where Y(t)={y(t ₁ ), y(t ₂ ), . . . , y(t _N )}={y ₁ , y ₂ , . . . , y _N }.

Example 1

This embodiment is a simulation experiment of the method of the present invention carried out under the MATLAB 2015a programming environment. The simulation algorithm of the magnetic resonance sounding signal extraction method based on frequency-selective singular spectrum analysis, referring to Figure 1, includes the following steps:

Step (1): Construct a pure MRS signal with a Larmor frequency of 2114Hz, an amplitude e ₀ of 180nV, a relaxation time T ₂ ^* of 0.1s, and a phase of 1.03. The signal sampling rate is 10kHz, the signal length is 500ms, and the data The number of points is 5000, as shown in Figure 3(a) and Figure 3(b). On the basis of this signal, 16 groups of noise-containing signals are simulated, and each group of signals is added with 100 power frequency interferences with an amplitude of 100nV, random phase, and an integer multiple of 50Hz between 0Hz and 5000Hz; the amplitude is 100nV, continuous Spike noise with a duration of 10ms and random noise with an amplitude of 60nV. The 16 groups of noises are statistically superimposed to form the observed MRS signal X ₁ (t) (row vector) with a signal-to-noise ratio of -13.3930dB. The time domain diagram is shown in Figure 3(a), and the power spectrum diagram is shown in Figure 3(b );

Step (2): pass the observed MRS signal X ₁ (t) through the bandpass filter shown in Figure 5, the adopted bandpass filter is a Chebyshev filter, the left boundary of the passband is 1900Hz, and the right boundary of the passband is The boundary is 2300Hz, the left boundary of the stopband is 1700Hz, the right boundary of the stopband is 2500Hz, the sideband area of the passband is attenuated by 0.1dB, and the cutoff area of the stopband is attenuated by 30dB, and the signal X ₂ (t) is obtained, as shown in Figure 4(a) .

Step (3): Calculate the power spectrum of the observed MRS signal X ₂ (t), as shown in Fig. 4(b). Make sure that the amplitude of the Larmor signal in the noisy signal ranks first;

Step (4): Embedding X ₂ (t) obtained through the broadband bandpass filter, the selected window length is Get the trajectory matrix H;

Decompose the trajectory matrix H by RSVD, let k=20, w=200, Satisfy the approximate conditions, get the singular values and singular vector U arranged in descending order, and draw the singular spectrum as shown in Figure 6;

According to the sorting position obtained in step (3), the Larmor signal amplitude is ranked first, so the first and second two singular values are selected for matrix reconstruction to obtain matrix X;

The matrix X is diagonally averaged to obtain the extracted MRS signal Y(t), the time domain diagram and power spectrum are shown in Figure 7(a) and Figure 7(b);

In order to verify the practicability of the method of the present invention, the signal-to-noise ratio (SNR) of the denoised MRS signal Y(t) is estimated. After calculation, its SNR=9.2147dB, which is 22.6076dB higher than the SNR before separation; then, envelope extraction and data fitting were performed on Y(t) to obtain the key parameters of the separation signal, the initial amplitude e ₀ and the relaxation time T ₂ ^* can be calculated, The relative errors are 1.6024% and 3.0142%, respectively, both of which are controlled within ±5%, meeting the application requirements.

Example 2

In this embodiment, the MRS signals collected in Shaoguo Town, Changchun City (the Larmor frequency is about 2332.5 Hz) are used as the processing object of the method of the present invention. As shown in Figure 1, the method for extracting magnetic resonance sounding signals based on frequency-selective singular spectrum analysis includes the following steps:

Step (1): A set of observed MRS signal X ₁ (t) (row vector) is collected by nuclear magnetic resonance sounding (MRS) water detector, as shown in Fig. 8(a) and Fig. 8(b), sampling The frequency is 25kHz, the number of data points is 12476, and the calculated signal-to-noise ratio is SNR=-7.1453dB;

Step (2): pass the observed MRS signal X ₁ (t) through the bandpass filter shown in Figure 9, the adopted bandpass filter is a Chebyshev filter, the left boundary of the passband is 2100Hz, and the right boundary of the passband The boundary is 2500Hz, the left boundary of the stopband is 1800Hz, the right boundary of the stopband is 2800Hz, the sideband area of the passband is attenuated by 0.1dB, and the cutoff area of the stopband is attenuated by 30dB, and the signal X ₂ (t) is obtained, as shown in Figure 8(a) .

Step (3): Calculate the power spectrum of the observed MRS signal X ₂ (t), as shown in blue in Figure 8(b). It is determined that the amplitude of the Larmor signal in the noisy signal is ranked second;

Step (4): Embedding X ₂ (t) obtained through the broadband bandpass filter, the selected window length is Get the trajectory matrix H;

Decompose the trajectory matrix H by RSVD, let k=20, w=200, Satisfy the approximate conditions, get the singular values and singular vector U arranged in descending order, and draw the singular spectrum as shown in Figure 10;

According to the sorting position obtained in step (3), the Larmor signal amplitude ranks second, so the third and fourth two singular values are selected for matrix reconstruction to obtain matrix X;

The matrix X is diagonally averaged to obtain the extracted MRS signal Y(t), the time domain diagram and power spectrum are shown in Figure 11(a) and Figure 11(b);

In order to verify the practicability of the method of the present invention, the signal-to-noise ratio (SNR) of the denoised MRS signal Y(t) is estimated. After calculation, its SNR=11.9542dB, which is 19.0995dB higher than the SNR before separation; then, envelope extraction and data fitting were performed on Y(t) to obtain the key parameters of the separation signal, the initial amplitude e ₀ and the relaxation time T ₂ ^* can be calculated, e ₀ = 63.7791nV, T ₂ ^* = 0.1890s, which is consistent with the actual hydrogeological drilling data.

Claims (5)

1. A magnetic resonance sounding signal extraction method based on frequency-selective singular spectrum analysis, is characterized in that, comprises the following steps: Step (1): Collect a set of observed MRS signals X ₁ (t) with a known Larmor frequency using the nuclear magnetic resonance sounding instrument; Step (2): pass the collected observed MRS signal X ₁ (t) through a broadband band-pass filter to obtain a signal X ₂ (t); Step (3): Perform power spectrum analysis on the signal X ₂ (t) obtained through the broadband bandpass filter, arrange the amplitudes of each frequency in descending order, and find the sorting position of the signal amplitude corresponding to the Larmor frequency; Step (4): performing singular spectrum analysis on the signal X ₂ (t) obtained through the broadband bandpass filter to extract the MRS signal, wherein: the singular spectrum analysis includes: Embedding to get the trajectory matrix H; Decompose the trajectory matrix H by RSVD to obtain singular values and singular vectors U arranged in descending order; According to the sorting position that step (3) obtains, select two singular values corresponding to it and carry out matrix reconstruction, obtain matrix C; The matrix C is diagonally averaged to obtain the extracted MRS signal Y(t). 2. according to a kind of magnetic resonance sounding signal extraction method based on frequency selection singular spectrum analysis according to claim 1, it is characterized in that, the concrete steps of embedding in the described step (4) are: The signal X ₂ (t) is a one-dimensional real sequence of length N: X ₂ (t)=(x ₁ ,x ₂ ,...,x _N ), the positive integer L is the length of the sliding window, 1<L<N , the value of L is determined by the following formula: P vectors are formed by embedding the original sequence signal X ₂ (t), and each vector is denoted by h _i h _i ＝(x _i ,x _i+1 ,…,x _i+L-1 ) ^T Where P=N-L+1, i=1,2,...,P, the result of the mapping forms a trajectory matrix H: 3. according to a kind of magnetic resonance sounding signal extracting method based on frequency selective singular spectrum analysis according to claim 1, it is characterized in that, the concrete steps of the RSVD decomposition in the described step (4) are: 1), set parameter k and parameter w, k is the k singular values taken for approximating the reconstruction matrix, w is a parameter used to ensure that the conditions of the reconstruction matrix are established, where w>k, w<L, w <P; 2), construct a Gaussian random matrix G _L×w with 0 mean and 1 variance; 3) Calculate the sampling matrix M _P×w of the trajectory matrix H _L×P : 4) Select the first k singular values of the sampling matrix M _P×w and perform SVD decomposition to obtain the left singular vector Q _P×k of the orthogonal matrix, the diagonal matrix Z _k×k , and the right singular vector Expressed as: 5), construct matrix T _L×k : _TL×k ＝H _L×P Q _P×k ; 6) Perform SVD decomposition on matrix T _L×k to obtain its left singular vector U _L×L , diagonal matrix Σ _L×k and right singular vector Expressed as: 7) Calculation matrix V _P×k : V _{P × k} = Q _{P × k} O _{k × k} ; 8), calculate the approximate reduced rank matrix of H _L×P Use RSVD decomposition to obtain the approximate reduced rank matrix of H, and the obtained approximate reduced rank matrix must meet the following conditions: λ _k+1 is the k+1th singular value of the trajectory matrix H; The diagonal matrix Σ _L×k is a matrix of order L×k, and its main diagonal elements are the k singular values of the MRS signal approximate reduced-rank matrix; the k singular values are arranged in descending order (λ ₁ ≥λ ₂ ≥... ≥λ _k ), the corresponding orthogonal singular vector set U′=(u ₁ ,u ₂ ,...u _k ), the contribution rate of each singular value is: Draw a singular spectrum according to the above formula. 4. according to a kind of magnetic resonance sounding signal extraction method based on frequency selective singular spectrum analysis according to claim 3, it is characterized in that, according to sorting position in described step (4), select two corresponding to it Singular values for matrix reconstruction include: According to step (3), find the sorting v of the magnitude of the MRS signal amplitude, and select the 2v-1 and 2v singular values in the singular spectrum; According to the selected singular value reconstruction matrix, the reconstruction steps are as follows: 1) Calculate the right singular vector W of the reconstructed signal matrix Among them, j=2v-1,2v; 2), reconstruct signal matrix C 5. according to a kind of magnetic resonance sounding signal extracting method based on frequency selective singular spectrum analysis according to claim 1, it is characterized in that, the concrete steps of the diagonal averaging in the described step (4) are: Singular spectrum analysis reconstructed signal matrix {c ₁ ,c ₂ ,…,c ^q ,…,c ^P } is transformed into corresponding reconstructed sequence {g ¹ ,g ² ,…,g ^q ,… by diagonal averaging ,g ^P }, where the sequence g ^q represents the qth singular spectrum analysis reconstruction sequence, The process is as follows: Among them, c ^q represents the qth singular spectrum analysis reconstructed signal matrix, is the element in row m and column n of matrix c ^q , Represents the dth element of the reconstructed sequence g ^q , q=1,2,...,P, d=1,2,...,N; the reconstructed sequence {g ₁ ,g ₂ ,...,g ^P } accumulate and sum to obtain the denoised MRS signal sequence Where Y(t)={y(t ₁ ), y(t ₂ ), . . . , y(t _N )}={y ₁ , y ₂ , . . . , y _N }.