阅读说明：本技术 一种磁共振波谱重建方法及系统 (Magnetic resonance spectrum reconstruction method and system ) 是由 郎俊 程达 于 2019-10-17 设计创作，主要内容包括：本发明涉及一种磁共振波谱重建方法及系统,该方法包括：首先获取欠采样矩阵X<Sub>0</Sub>；然后对欠采样矩阵X<Sub>0</Sub>进行奇异值分解,得到X<Sub>0</Sub>＝UV<Sup>T</Sup>；再将U和V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型；接着引入中间变量将磁共振波谱低秩模型转化为最小秩优化问题；采用交替方向乘子法求解最小秩优化问题,得到U<Sup>k+1</Sup>和V<Sup>k+1</Sup>,并补全矩阵X：最后当X<Sup>k+1</Sup>为二维磁共振信号时,对X<Sup>k+1</Sup>做二维傅里叶变换得到完整的二维磁共振波谱。本发明利用磁共振信号可以进行范德蒙德分解的性质,对欠采样矩阵X<Sub>0</Sub>进行奇异值分解得矩阵U和V,通过采用交替方向乘子法分别求解U和V,完成对U和V的重建,从而实现对欠采样矩阵X<Sub>0</Sub>的补全,重建速度快且精度高。(The invention relates to a magnetic resonance spectrum reconstruction method and a system, wherein the method comprises the following steps: firstly, an undersampling matrix X is obtained 0 (ii) a Then, the undersampling matrix X is processed 0 Singular value decomposition is carried out to obtain X 0 ＝UV T (ii) a Converting each column of U and V into a Hankel matrix, and constructing a magnetic resonance spectrum low-rank model; then, introducing an intermediate variable to convert the magnetic resonance spectrum low-rank model into a minimum-rank optimization problem; solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain U k+1 And V k+1 And completing the matrix X: finally when X is k+1 For two-dimensional magnetic resonance signals, for X k+1 And performing two-dimensional Fourier transform to obtain a complete two-dimensional magnetic resonance spectrum. The invention utilizes the characteristic that the magnetic resonance signal can be subjected to Van der Mond decomposition to the undersampled matrix X 0 Singular value decomposition is carried out to obtain matrixes U and V, the U and the V are respectively solved by adopting an alternative direction multiplier method, the U and the V are reconstructed, and therefore the undersampled matrix X is realized 0 The completion, the reconstruction speed is fast and the precision is high.)

1. A magnetic resonance spectroscopy reconstruction method, characterized in that the magnetic resonance spectroscopy reconstruction method specifically comprises the steps of:

s1: obtaining an undersampling matrix X₀；

S2: for the undersampling matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Is a right singular vector, sigma is a singular value;

taking the undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^TThen X₀＝UV^T；

S3: converting each column of the U and the V into a Hankel matrix, and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is an operation operator for converting the one-dimensional data into a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;

s4: introducing intermediate variables Bi ═ RQIU and Ci ═ RQIV, and converting the magnetic resonance spectrum low-rank model into a minimum rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain the U^k+1And said V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1Obtaining a matrix U after k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

s5: when said X is^k+1For the two-dimensional magnetic resonance signal, for the X^k+1And performing two-dimensional Fourier transform to obtain a complete two-dimensional magnetic resonance spectrum.

2. The method according to claim 1, wherein in step S2, r is within 6 times the number of peaks in the actual two-dimensional magnetic resonance spectrum.

3. The magnetic resonance spectroscopy reconstruction method as set forth in claim 1, wherein in step S4, specifically:

s41: converting the minimum rank optimization problem into an augmented Lagrangian function:

s.t.P_Ω(X)＝X₀

wherein Di, Mi is Lagrange multiplier, mu is regularization parameter, mu is 0, <, > represents the inner product of the matrix in Hilbert space, < A, B > is trace (AB), trace (AB) represents the trace of the matrix;

s42: and alternately solving the following subproblems by adopting an alternate direction multiplier method for the augmented Lagrange function to obtain:

s.t.P_Ω(X)＝X₀

s43: u is obtained by solving the following least squares problem:

s44: solving the following least squares problem yields V:

s45: iteratively updating X by:

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

iteratively updating B by the following formula_i、C_i：

Wherein k is the number of iterations; r is an accompanying operator of the operation operator R; qi is the accompanying operator of Qi;

4. A method of magnetic resonance spectroscopy as claimed in claim 3, wherein the singular value reduction operator

said C is_iThe optimization sub-problem of (a) is expressed as:

wherein σ_BiRepresenting the singular values of the matrix Bi; sigma_CiSingular values representing the matrix Ci; i | · | purple wind₀Representing the 0 norm of the matrix.

5. A method of magnetic resonance spectroscopy as claimed in claim 3 wherein the value of μ increases progressively with each iteration, and when an iteration stop criterion is reached, the value of β is increased and the iteration is continued until β reaches a set value and the iteration ends.

6. The method according to claim 5, wherein the initial β value is set to 20-30, the initial μ value is set to 0.005-0.02, when one iteration is completed, the value of μ in the next iteration is 1.02-1.12 times the value of μ in the previous iteration, when the iteration stop criterion is met, the value of β in the next iteration is twice the value of β in the previous iteration, and the iteration is continued until β reaches a maximum value of 2²²-2⁴²Ending the iteration; the iteration stop criterion is set to the matrix X completed in two adjacent iterations_k+1And X_kIs less than a threshold η, η is set to 10^-7-10^-5。

7. The magnetic resonance spectroscopy reconstruction method as set forth in claim 3, wherein in step S43, specifically:

R^*the companion operator, which is the operator R, is defined as follows:

R^*d converts a matrix D into a vector,

qi is the companion operator of Qi, defined as follows:

d is a vector and d is a vector,

therefore, the first and second electrodes are formed on the substrate,

wherein C represents a complex field, and M, N represents the number of rows and columns of the matrix respectively; w1 is an mx 1 vector and its k-th element is the number of elements on the k-th sub-diagonal of said matrix RQiU, w2 is an nx1 vector and its k-th element is the number of elements on the k-th sub-diagonal of said matrix RQiV; defining a matrix T1 to be r rows w1, and T2 to be r rows w 2;

order to

Therefore, the first and second electrodes are formed on the substrate,

conversion to:

and iterate U by updating each line of U by the following formula:

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

wherein, U_(a,:) ^k+1Representing the matrix U generated after k +1 iterations of the matrix U^k+1The elements of line a of (a); y1_(a,:)Is an element of row a of matrix Y1; t is_aAs a diagonal matrix, T_aThe main diagonal element is T1_(a,:)，T1_(a,:)Is an element of row a of the matrix T1.

8. The magnetic resonance spectroscopy reconstruction method as set forth in claim 7, wherein the step S44 is specifically:

order to

Obtaining:

is converted into

And iterating V by updating each line of V by the following equation:

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

wherein, V_(b,:) ^k+1Representing the matrix V generated after k +1 iterations of the matrix V^K+1The element of the nth row; y2_(b,:)Is an element of the nth row of matrix Y2; t is_bAs a diagonal matrix, T_bThe main diagonal element is T2_(b,:)，T2_(b,:)Is an element of the b-th row of the matrix T2.

9. The magnetic resonance spectroscopy reconstruction method as set forth in claim 1, wherein the step S5 further includes:

and when the three-dimensional magnetic resonance signal is obtained, respectively carrying out two-dimensional Fourier transform on all two-dimensional magnetic resonance signals forming the three-dimensional magnetic resonance signal, and carrying out three-dimensional Fourier transform on the three-dimensional magnetic resonance signal to obtain a complete three-dimensional magnetic resonance spectrum.

10. A magnetic resonance spectroscopy reconstruction system, the two-dimensional magnetic resonance spectroscopy reconstruction system comprising: the system comprises an acquisition module, a singular value decomposition module, a low-rank model module, an optimization module and a transformation module;

the acquisition module is used for acquiring an undersampling matrix X₀；

The singular value decomposition module is used for undersampling the matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Is a right singular vector, sigma is a singular value;

taking the undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, and let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^TThen X₀＝UV^T；

The low-rank model module is used for converting each column of the U and the V into a Hankel matrix and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is an operation operator for converting the one-dimensional data into a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;

the optimization module is used for introducing intermediate variables Bi ═ RQIU and Ci ═ RQIV and converting the magnetic resonance spectrum low-rank model into a minimum-rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain the U^k+1And said V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1Obtaining a matrix U after k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

the transformation module is used for converting the X^k+1For the two-dimensional magnetic resonance signal, for the X^k+1And performing two-dimensional Fourier transform to obtain a two-dimensional frequency spectrum of the spectrum.

Technical Field

The invention relates to the field of magnetic resonance spectrum reconstruction, in particular to a magnetic resonance spectrum reconstruction method and a magnetic resonance spectrum reconstruction system.

Background

Nmr spectroscopy, which is one of the most powerful tools for the qualitative and sometimes quantitative analysis of the composition and structure of various organic and inorganic substances, is the study of the absorption of radio frequency radiation by atomic nuclei. The nuclear magnetic resonance technology can provide information of the chemical structure and molecular dynamics of molecules, becomes a conventional technical means for molecular structure analysis and physical and chemical property characterization, can determine the environment of almost all common functional groups, has strong intuition, particularly can directly reflect the molecular skeleton through a carbon spectrum, is easy to explain a spectrogram, is widely applied to the fields of physics, chemistry, biology, medicine, food and the like, and is an indispensable means for conventional analysis in chemistry.

In practical application, due to higher experimental cost and hardware limitation and other inevitable reasons, it is difficult to completely sample a two-dimensional or three-dimensional magnetic resonance free attenuation signal, so that it is important to undersample a magnetic resonance signal and reconstruct a complete magnetic resonance signal from data obtained by undersampling. Different reconstruction methods have great differences in reconstruction accuracy, minimum undersampling rate and reconstruction speed, while the existing methods all require long reconstruction time, and particularly the time required for reconstructing a three-dimensional magnetic resonance spectrum is increased in proportion. Therefore, it is a problem to be solved currently to provide a new reconstruction method with faster reconstruction speed and higher reconstruction accuracy.

Disclosure of Invention

Technical problem to be solved

The invention provides a magnetic resonance spectrum reconstruction method, and aims to provide a new reconstruction method with higher reconstruction speed and higher reconstruction precision.

(II) technical scheme

In order to solve the above problems, the present invention provides a magnetic resonance spectroscopy reconstruction method, including the steps of:

s1: obtaining an undersampling matrix X₀；

S2: for the undersampling matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Is a right singular vector, sigma is a singular value;

taking the undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^TThen X₀＝UV^T；

S3: converting each column of the U and the V into a Hankel matrix, and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is an operation operator for converting the one-dimensional data into a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;the square of the F-norm representing the matrix; beta is a regularization parameter used for balancing the importance of the two terms, beta>0; x represents the completed matrix; the set Ω represents a set of positions of sampling points;

i.e. P_ΩHadamard product representation with X as template P_ΩUndersampling the X, denoted as P_Ω(X), the matrix X corresponds to P_ΩThe element of the position with the number 1 remains unchanged, and the element of the position corresponding to the template number 0 is set to 0;

s4: introducing intermediate variables Bi ═ RQIU and Ci ═ RQIV, and converting the magnetic resonance spectrum low-rank model into a minimum rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain the U^k+1And said V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1Obtaining a matrix U after k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

s5: when said X is^k+1For the two-dimensional magnetic resonance signal, for the X^k+1And performing two-dimensional Fourier transform to obtain a complete two-dimensional magnetic resonance spectrum.

Preferably, in the step S2, the r is within 6 times of the number of peaks of the actual two-dimensional magnetic resonance spectrum.

Preferably, in step S4, specifically, the method includes:

s41: converting the minimum rank optimization problem into an augmented Lagrangian function:

s.t.P_Ω(X)＝X₀

wherein Di, Mi is Lagrange multiplier, mu is regularization parameter, mu is 0, <, > represents the inner product of the matrix in Hilbert space, < A, B > is trace (AB), trace (AB) represents the trace of the matrix;

s42: and alternately solving the following subproblems by adopting an alternate direction multiplier method for the augmented Lagrange function to obtain:

s.t.P_Ω(X)＝X₀

s43: u is obtained by solving the following least squares problem:

s44: solving the following least squares problem yields V:

s45: iteratively updating X by:

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

iteratively updating B by the following formula_i、C_i：

Wherein k is the number of iterations; r is an accompanying operator of the operation operator R; qi isA companion operator of Qi;

the operator is reduced for the singular value and,

is the said B_iAs a result of k +1 iterations,

is the C_iResult after k +1 iterations, μ^kIs the result of k iterations of said mu.

Preferably, the singular value reduction operator

For the hard threshold operator, the rank of the matrix is 0 norm of the singular value, in said step S42, B_iThe optimization sub-problem of (a) is expressed as:

said C is_iThe optimization sub-problem of (a) is expressed as:

wherein σ_BiRepresenting the singular values of the matrix Bi; sigma_CiSingular values representing the matrix Ci; i | · | purple wind₀Representing the 0 norm of the matrix.

Preferably, each iteration, the value of μ is gradually increased, and when an iteration stop criterion is reached, the value of β is increased and the iteration is continued until the value of β reaches a set value, and the iteration is ended.

Preferably, the initial value of β is set to 20-30, the initial value of μ is 0.005-0.02, after one iteration is completed, the value of μ in the next iteration is 1.02-1.12 times the value of μ in the last iteration, and when the iteration stop criterion is reached, the value of β in the next iteration is the value of β in the last iterationUntil said beta reaches a maximum value of 2²²-2⁴²Ending the iteration; the iteration stop criterion is set to the matrix X completed in two adjacent iterations_k+1And X_kIs less than a threshold η, η is set to 10^-7-10^-5。

Preferably, in step S43, specifically, the method includes:

R^*the companion operator, which is the operator R, is defined as follows:

R^*d converts a matrix D into a vector,

the first element of the vector is the sum of the elements on the first sublevel of the matrix D;

qi is the companion operator of Qi, defined as follows:

d is a vector and d is a vector,indicating that the elements of vector d are placed in the matrix Qi^*d]In the p columns of (1);

therefore, the first and second electrodes are formed on the substrate,

wherein C represents a complex field, and M, N represents the number of rows and columns of the matrix respectively; w1 is an mx 1 vector and its k-th element is the number of elements on the k-th sub-diagonal of said matrix RQiU, w2 is an nx1 vector and its k-th element is the number of elements on the k-th sub-diagonal of said matrix RQiV; defining a matrix T1 to be r rows w1, and T2 to be r rows w 2;

order to

Therefore, the first and second electrodes are formed on the substrate,

conversion to:

and iterate U by updating each line of U by the following formula:

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

wherein, U_(a,:) ^k+1Representing the matrix U generated after k +1 iterations of the matrix U^k+1The elements of line a of (a); y1_(a,:)Is an element of row a of matrix Y1; t is_aAs a diagonal matrix, T_aThe main diagonal element is T1_(a,:)，T1_(a,:)Is an element of row a of the matrix T1.

Preferably, the step S44 is specifically:

order to

Obtaining:

is converted into

And iterating V by updating each line of V by the following equation:

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

wherein, V_(b,:) ^k+1Representing the matrix V generated after k +1 iterations of the matrix V^K+1The element of the nth row; y2_(b,:)Is an element of the nth row of matrix Y2; t is_bAs a diagonal matrix, T_bThe main diagonal element is T2_(b,:)，T2_(b,:)Is an element of the b-th row of the matrix T2.

Preferably, the step S5 further includes:

and when the three-dimensional magnetic resonance signal is obtained, respectively carrying out two-dimensional Fourier transform on all two-dimensional magnetic resonance signals forming the three-dimensional magnetic resonance signal, and carrying out three-dimensional Fourier transform on the three-dimensional magnetic resonance signal to obtain a complete three-dimensional magnetic resonance spectrum.

Preferably, the present invention also provides a magnetic resonance spectroscopy reconstruction system comprising: the system comprises an acquisition module, a singular value decomposition module, a low-rank model module, an optimization module and a transformation module;

the acquisition module is used for acquiring an undersampling matrix X₀；

The singular value decomposition module is used for undersampling the matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Is a right singular vector, sigma is a singular value;

taking the undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, and let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^TThen X₀＝UV^T；

The low-rank model module is used for converting each column of the U and the V into a Hankel matrix and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is one-dimensional dataConverting into an operation operator of a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;the square of the F-norm representing the matrix; beta is a regularization parameter used for balancing the importance of the two terms, beta>0; x represents the completed matrix; the set Ω represents a set of positions of sampling points;

i.e. P_ΩHadamard product representation with X as template P_ΩUndersampling said X, denoted as P Ω (X), such that said matrix X corresponds to P_ΩThe element of the position with the number 1 remains unchanged, and the element of the position corresponding to the template number 0 is set to 0;

the optimization module is used for introducing intermediate variables Bi ═ RQIU and Ci ═ RQIV and converting the magnetic resonance spectrum low-rank model into a minimum-rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain the U^k+1And said V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1Obtaining a matrix U after k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

the transformation module is used for converting the X^k+1For the two-dimensional magnetic resonance signal, for the X^k+1And performing two-dimensional Fourier transform to obtain a two-dimensional frequency spectrum of the spectrum.

(III) advantageous effects

The method comprises the steps of carrying out singular value decomposition on an undersampled matrix, converting each column of the matrix obtained by decomposition into a Hankel matrix, establishing a magnetic resonance spectrum low-rank model, obtaining a reconstructed matrix by adopting an alternative multiplier method, and finally carrying out Fourier transform on the reconstructed matrix to obtain a complete magnetic resonance spectrum. The reconstruction method has the advantages of higher reconstruction speed and higher reconstruction precision.

Drawings

FIG. 1 is a flow chart of a magnetic resonance spectroscopy reconstruction method of the present invention;

FIG. 2 is a flow chart of solving a minimum rank optimization problem in the present invention;

FIG. 3 is a schematic diagram of a magnetic resonance spectroscopy reconstruction system of the present invention;

FIG. 4 is a full sample two-dimensional magnetic resonance spectrum of 256 x 128 in size according to an embodiment of the present invention;

FIG. 5 is a magnetic resonance spectrum reconstructed using the method of the present invention to complement an undersampled matrix;

FIG. 6 is a comparison of the present invention with other methods for reconstructing results at 20% sampling rate;

FIG. 7 is a comparison of the results of reconstruction at 20% sampling rate using different threshold reduction operators.

[ description of reference ]

1: an acquisition module; 2: a singular value decomposition module; 3: a low rank model module; 4: an optimization module; 5: and a transformation module.

Detailed Description

For the purpose of better explaining the present invention and to facilitate understanding, the present invention will be described in detail by way of specific embodiments with reference to the accompanying drawings.

It should be noted that all the directional indicators (such as up, down, left, right, front, and rear … …) in the embodiment of the present invention are only used to explain the relative position relationship between the components, the movement situation, etc. in a specific posture (as shown in the drawing), and if the specific posture is changed, the directional indicator is changed accordingly.

In addition, the descriptions related to "first", "second", etc. in the present invention are only for descriptive purposes and are not to be construed as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one such feature. In the description of the present invention, "a plurality" means at least two, e.g., two, three, etc., unless specifically limited otherwise.

In the present invention, unless otherwise expressly stated or limited, the terms "connected," "secured," and the like are to be construed broadly, and for example, "secured" may be a fixed connection, a removable connection, or an integral part; can be mechanically or electrically connected; they may be directly connected or indirectly connected through intervening media, or they may be connected internally or in any other suitable relationship, unless expressly stated otherwise. The specific meanings of the above terms in the present invention can be understood by those skilled in the art according to specific situations.

The invention provides a magnetic resonance spectrum reconstruction method, which is shown in figure 1: the flow chart of the magnetic resonance spectrum reconstruction method of the invention is shown, and the method specifically comprises the following steps:

s1: obtaining an undersampling matrix X₀。X₀Is an M × N matrix, M, N respectively represents the matrix X₀The number of rows and columns.

S2: to the undersampling matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Right singular vectors, Σ singular values.

Taking undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^TThen X₀＝UV^T. U is an M × r matrix and V is an N × r matrix.

The one-dimensional magnetic resonance signal can be represented as a superposition of several complex exponential decay signals, while the two-dimensional magnetic resonance signal is a superposition of several complex exponential decay signals multiplied by two, and the m-th row and n-th column elements in a matrix P formed by the two-dimensional magnetic resonance signal are represented as:wherein, y_i＝exp(j2πf_1i-τ_1i),z_i＝exp(j2πf_2i-τ_2i) F represents normalized frequency, tau is attenuation factor, b represents that the two-dimensional magnetic resonance spectrum has g spectrum peaks; a complete two-dimensional magnetic resonance signal matrix P can be decomposed into two van der mond matrices as follows:

D＝diag[d₁,d₂，…d_g]

P＝YDZ^Tdividing P into D ═ D_LD_RTwo diagonal matrices, multiplied by Y and Z respectively^TIs provided with U₁＝YD_L,V₁ ^T＝D_RZ^TThen X is equal to U₁V₁ ^T,U₁、V₁Is a van der mond matrix.

Undersampling matrix X₀Also obtained by undersampling the magnetic resonance signals, in the case of matrix X₀The singular value decomposition yields U, V if we enable U, V to approximate the Van der Mond matrix U₁、V₁Then the complete magnetic resonance signal can be reconstructed, and the application utilizes the property that the magnetic resonance signal can be decomposed into two Van der Monde matrixes, and the complementary matrix X is obtained₀Turn to U, V approaching U₁、V₁Therefore, the completion method is more accurate and faster.

S3: converting each column of U and V into a Hankel matrix, and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is an operation operator for converting the one-dimensional data into a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;

the square of the F-norm representing the matrix; beta is a regularization parameter used for balancing the importance of the two terms, beta>0; x represents the completed matrix; the set Ω represents a set of positions of sampling points;

i.e. P_ΩHadamard product representation with X as template P_ΩX is undersampled and denoted as P Ω (X), and matrix X corresponds to P_ΩThe element of the position with the number 1 remains unchanged, and the element of the position corresponding to the template number 0 is set to 0;

assuming e is data of dimension 1 × a, Re is:

where k is set to 0.1A, and e (f) represents the f-th element in e.

S4: introducing intermediate variables Bi ═ RQIU and Ci ═ RQIV, and converting a magnetic resonance spectrum low-rank model into a minimum-rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain U^k+1And V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1The matrix U is subjected to k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

s5: when X is present^k+1For two-dimensional magnetic resonance signals, for X^k+1And performing two-dimensional Fourier transform to obtain a complete two-dimensional magnetic resonance spectrum.

Further, in step S2, r is within 6 times the number of peaks of the actual two-dimensional magnetic resonance spectrum. In some practical situations, the number of spectral peaks of the magnetic resonance signal is not determined, so that the spectral peak r needs to be preset, and the preset number can be successfully reconstructed within 6 times of the actual number of spectral peaks, namely r is less than 6 g.

As shown in fig. 2: as shown in the flowchart of solving the minimum rank optimization problem in the present invention, in a preferred embodiment, step S4 specifically includes:

s41: converting the minimum rank optimization problem into an augmented Lagrangian function:

s.t.P_Ω(X)＝X₀

wherein Di, Mi is Lagrange multiplier, mu is regularization parameter, mu >0, < - > represents the inner product of the matrix in Hilbert space, < G, H > -trace (GH), trace (·) represents the trace of the matrix;

s42: and (3) alternately solving the following subproblems by adopting an alternate direction multiplier method for the augmented Lagrange function to obtain:

s.t.P_Ω(X)＝X₀

s43: u is obtained by solving the following least squares problem:

s44: solving the following least squares problem yields V:

s45: iteratively updating X by:

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

iteratively updating B by the following formula_i、C_i：

Wherein k is the number of iterations; r is an accompanying operator of the operation operator R; qi is the accompanying operator of Qi;

the operator is reduced for the singular value and,

is B_iAs a result of k +1 iterations,

is C_iResult after k +1 iterations, μ^kIs the result of k iterations of mu.

In a more preferred embodiment, the singular value reduction operator

For the hard threshold operator, the rank of the matrix is 0 norm of the singular value, in step S42, B_iThe optimization sub-problem of (a) is expressed as:

C_ithe optimization sub-problem of (a) is expressed as:

wherein σ_BiRepresenting the singular values of the matrix Bi; sigma_CiSingular values representing the matrix Ci; i | · | purple wind₀Representing the 0 norm of the matrix.

Each RQIU and RQIV is a Hankel matrix and is a rank-1 matrix, most information is in the first singular value for rank-1, main information can be better reserved by adopting iteration of a hard threshold operator, interference is removed, the reconstruction speed and the accuracy effect are better under the condition of low sampling rate, and the hard threshold operator is more suitable for the reconstruction problem. The matrix adopts a 0 norm model and uses a hard threshold operator for iteration, and the effect is better than that of other threshold reduction operators.

And when the iteration stopping criterion is reached, increasing beta and continuing the iteration until the beta reaches a set value, and finishing the iteration. Setting the initial value of beta to be 20-30, the initial value of mu to be 0.005-0.02, when one iteration is finished, the value of mu in the next iteration is 1.02-1.12 times of the value of mu in the last iteration, when the iteration stop criterion is reached, the value of beta in the next iteration is twice of the value of beta in the last iteration, and continuing the iteration until the beta reaches the maximum value 2²²-2⁴²Ending the iteration; the iteration stop criterion is set to the matrix X completed in two adjacent iterations_k+1And X_kIs less than a threshold η, η is set to 10^-7-10^-5。

Further, step S43 is specifically:

R^*the companion operator, which is the operator R, is defined as follows:

R^*d converts a matrix D into a vector,

the first element of the vector is the sum of the elements on the first sublevel of the matrix D;

qi is the companion operator of Qi, defined as follows:

d is a vector and d is a vector,

indicating that the elements of vector d are placed in the matrix Qi^*d]In the p columns of (1);

therefore, the first and second electrodes are formed on the substrate,

wherein C represents a complex field, and M, N represents the number of rows and columns of the matrix respectively; w1 is an mx 1 vector and its k-th element is the number of elements on the kth sub-diagonal of the matrix RQiU, w2 is an nx1 vector and its k-th element is the number of elements on the kth sub-diagonal of the matrix RQiV; defining a matrix T1 to be r rows w1, and T2 to be r rows w 2;

order to

Therefore, the first and second electrodes are formed on the substrate,

conversion to:

and iterate U by updating each line of U by the following formula:

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

wherein, U_(a,:) ^k+1Representing the matrix U generated after k +1 iterations of the matrix U^k+1The elements of line a of (a); y1_(a,:)Is an element of row a of matrix Y1; t is_aAs a diagonal matrix, T_aThe main diagonal element is T1_(a,:)，T1_(a,:)Is an element of row a of the matrix T1.

Further, step S44 is specifically:

order to

Obtaining:

is converted into

And iterating V by updating each line of V by the following equation:

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

wherein, V_(b,:) ^k+1Representing the matrix V generated after k +1 iterations of the matrix V^K+1The element of the nth row; y2_(b,:)Is an element of the nth row of matrix Y2; t is_bAs a diagonal matrix, T_bThe main diagonal element is T2_(b,:)，T2_(b,:)Is an element of the b-th row of the matrix T2.

Finally, step S5 further includes:

when the three-dimensional magnetic resonance signal is obtained, two-dimensional Fourier transform is respectively carried out on all two-dimensional magnetic resonance signals forming the three-dimensional magnetic resonance signal, and three-dimensional Fourier transform is carried out on the three-dimensional magnetic resonance signal to obtain a complete three-dimensional magnetic resonance spectrum.

The invention also provides a magnetic resonance spectrum reconstruction system, as shown in fig. 3: a schematic diagram of a magnetic resonance spectroscopy reconstruction system of the present invention is shown, including: the system comprises an acquisition module 1, a singular value decomposition module 2, a low-rank model module 3, an optimization module 4 and a transformation module 5;

the acquisition module 1 is used for acquiring an undersampled matrix X₀；

Singular value decomposition module 2 for undersampling matrix X₀Singular value decomposition is carried out to obtain X₀＝U₀ΣV₀ ^TWherein, U₀Is the left singular vector, V₀Is a right singular vector, sigma is a singular value;

taking undersampling matrix X₀The first r singular values and the left singular vectors and the right singular vectors corresponding to the r singular values, and let Σ ═ Σ_LΣ_R，U＝U₀Σ_L，VT＝Σ_RV₀ ^TThen X₀＝UVT；

The low-rank model module 3 is used for converting each column of U and V into a hankel matrix and constructing a magnetic resonance spectrum low-rank model:

s.t.P_Ω(X)＝X₀

where r represents the number of columns of matrices U and V; r is an operation operator for converting the one-dimensional data into a Hankel matrix; the operator Qi is used for extracting the ith column of the matrix; rank (·) denotes the rank of the matrix;

the square of the F-norm representing the matrix; beta is a regularization parameter used for balancing the importance of the two terms, beta>0; x represents the completed matrix; the set Ω represents a set of positions of sampling points;

i.e. P_ΩHadamard product representation with X as template P_ΩX is undersampled, denoted P Ω (X), such that matrix X corresponds to P_ΩThe element of the position with the number 1 remains unchanged, and the element of the position corresponding to the template number 0 is set to 0;

the optimization module 4 is configured to introduce intermediate variables Bi ═ RQiU and Ci ═ RQiV, and convert the magnetic resonance spectrum low-rank model into a minimum rank optimization problem:

s.t.Bi＝RQiU,Ci＝RQiV,P_Ω(X)＝X₀

solving the minimum rank optimization problem by adopting an alternating direction multiplier method to obtain U^k+1And V^k+1And completing the matrix X:

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

wherein X^k+1The matrix X is subjected to k +1 times of iteration; u shape^k+1The matrix U is subjected to k +1 iterations; v^k+1The matrix V is subjected to k +1 iterations; p_ΩcIs P_ΩComplement of (2), its value is equal to P_ΩIn contrast, i.e. P_ΩIs a position of 1, P_Ωc0 at the corresponding position; p_ΩA position of 0, P_ΩcAt the corresponding position is 1;

the conversion module 5 is used for X_k+1For two-dimensional magnetic resonance signals, for X_k+1Two-dimensional Fourier transform is carried out to obtain a two-dimensional frequency spectrum

Now, the undersampled matrix X is illustrated by the specific example₀The size of the complete two-dimensional magnetic resonance spectroscopy data is 256 × 128, and contains 15 spectral peaks, and the two-dimensional magnetic resonance spectroscopy is shown in fig. 4. The reconstruction method of the magnetic resonance spectrum provided by the invention reconstructs X₀And obtaining reconstructed data X, and performing two-dimensional fourier transform on X to obtain a frequency spectrum thereof, as shown in fig. 5.