A magnetic resonance spectroscopy reconstruction method based on simultaneous row and column low-rank constraints

Abstract

The present invention discloses a magnetic resonance spectrum reconstruction method based on simultaneous low-rank constraints on rows and columns. It belongs to the field of digital image processing technology. It is a magnetic resonance spectrum reconstruction method that uses the binary exponential structure of a two-dimensional magnetic resonance spectrum to perform low-rank constraints on the Hankel matrices composed of its rows and columns. First, each row and each column of the two-dimensional magnetic resonance spectrum is constructed into a corresponding Hankel matrix, and then a non-convex norm is used to constrain the low-rank characteristics of the Hankel matrix, and finally the magnetic resonance spectrum reconstruction model established is optimized and solved by the alternating direction multiplier method. The present invention applies a non-convex low-rank constraint to the two-dimensional magnetic resonance spectrum, effectively extracts its binary exponential structure, and simultaneously accurately solves each sub-problem of the reconstruction model, so that the magnetic resonance spectrum reconstructed from the undersampled spatiotemporal encoding data effectively suppresses the artifact phenomenon, and the spectral peak in the reconstructed spectrum is closer to the real spectral peak, so it can be used for the reconstruction of the magnetic resonance spectrum.

Description

Technical Field

[0001] The present invention belongs to the technical field of digital image processing, and particularly relates to constructing a Hankel matrix for a magnetic resonance spectrum and performing non-convex low-rank constraints to achieve reconstruction, which is used for high-resolution restoration of the magnetic resonance spectrum. Background Art

[0002] Magnetic resonance spectroscopy is one of the important tools for analyzing the composition, structure and interaction of substances. It is widely used in chemical and biological engineering, but its high dimensionality and resolution lead to a sharp increase in data acquisition time. Spatiotemporal coding technology can encode the frequency information of magnetic resonance spectroscopy in the indirect dimension in parallel, thereby speeding up data acquisition. However, in order to reduce the requirements for hardware equipment and noise levels, it is necessary to implement non-uniform undersampling of spatiotemporal coded data, resulting in the loss of some data. On this basis, establishing an appropriate magnetic resonance spectroscopy reconstruction model and designing a feasible optimization algorithm are the key to reconstructing high-resolution magnetic resonance spectroscopy from undersampled spatiotemporal coded data.

[0003] Since the two-dimensional magnetic resonance spectrum has a binary exponential structure, constructing a Hankel matrix for the magnetic resonance spectrum and constraining its low-rank characteristics is more in line with the signal prior information than constraining the sparsity of the magnetic resonance spectrum, thereby improving the quality of the reconstructed magnetic resonance spectrum. However, the existing reconstruction methods generally only construct the Hankel matrix for the rows or columns of the two-dimensional magnetic resonance spectrum, which makes it difficult to fully utilize the characteristics of the signal itself. At the same time, the commonly used nuclear norm cannot impose strong low-rank constraints on the Hankel matrix, which ultimately limits the reconstruction performance of the magnetic resonance spectrum.

[0004] In fact, constructing the Hankel matrix for the rows and columns of the two-dimensional magnetic resonance spectrum and constraining its low rank property can fully utilize the binary exponential structure. ε The norm is closer to the rank of the matrix than the nuclear norm and can effectively constrain the low-rank characteristics of the matrix. Therefore, h can be used ε The norm simultaneously constrains the Hankel matrix constructed by the rows and columns of the magnetic resonance spectrum, and accurately solves the corresponding reconstruction model, making the reconstructed magnetic resonance spectrum closer to the real magnetic resonance spectrum. Summary of the invention

[0005] The purpose of the present invention is to make full use of the binary index structure of the two-dimensional magnetic resonance spectrum and propose a magnetic resonance spectrum reconstruction method based on simultaneous low-rank constraints on rows and columns. The method constructs Hankel matrices for the rows and columns of the two-dimensional magnetic resonance spectrum respectively, and uses a non-convex norm to constrain its low-rank characteristics, and then accurately solves the sub-problems about each variable in the reconstruction model under the optimization framework of the alternating direction multiplier method, so that the obtained reconstructed magnetic resonance spectrum effectively suppresses artifacts, and the reconstructed spectrum peak is closer to the real spectrum peak.

[0006] The specific steps include:

[0007] (1) Input a two-dimensional magnetic resonance spectroscopy signal's spatiotemporal encoding data y and an undersampling matrix U. For a two-dimensional magnetic resonance spectroscopy signal X to be reconstructed with N row vectors and M column vectors, the element in the mth column of the nth row is represented as x n，m ;

[0008] (2) For the vector x in the mth column of X m , construct its corresponding Hankel matrix It can be expressed as:

[0009]

[0010] in, represents the linear operator that transforms a column vector into a Hankel matrix, and L represents the size parameter;

[0011] (3) For the vector x in the nth row of X n , construct its corresponding Hankel matrix It can be expressed as:

[0012]

[0013] in, represents the linear operator that transforms a row vector into a Hankel matrix, and K represents the size parameter;

[0014] (4) Using non-convex h ε The low-rank property of the Hankel matrix constructed with norm constraints, h represents any matrix L ε The norm is defined as:

[0015]

[0016] Among them, σ r is the rth singular value of L, ε>0 is a small constant;

[0017] (5) When using h ε Based on the low-rank characteristics of the Hankel matrix constructed by simultaneously constraining the rows and columns in X, a reconstruction model for two-dimensional magnetic resonance spectroscopy is established:

[0018]

[0019] Among them, F M represents the Fourier transform matrix of size M×M, I N represents the identity matrix of size N×N, represents the Kronecker product operator of matrices, vec(·) represents the operator that stacks the columns of a matrix into a vector, represents the square of the second norm of the vector, λ represents the regularization parameter of the spatiotemporal coding data fidelity term;

[0020] (6) Introducing intermediate variables into the reconstruction model and And the sub-problems of each optimization variable in the reconstruction model are solved according to the alternating direction multiplier method:

[0021] (6a) Fix X and R n , find the solution about L m Sub-questions:

[0022]

[0023] Among them, β represents the penalty parameter, T m Indicates for L m The Lagrange multiplier introduced, Represents the square of the Frobenius norm of the matrix;

[0024] (6b) Fix X and L m , find the solution about R n Sub-questions:

[0025]

[0026] Among them, Z n Indicates that for R n The Lagrange multiplier introduced;

[0027] (6c) Fixed L m and R n , solve the subproblem about X:

[0028]

[0029] (7) After obtaining the solution of each optimization variable, update the Lagrange multiplier T m and Z n , and the penalty parameter β:

[0030]

[0031]

[0032] β=τβ

[0033] Where τ>1 represents a growth factor, and steps (6) to (7) are repeated until the reconstructed magnetic resonance spectrum meets the conditions or the number of iterations reaches a preset upper limit.

[0034] The innovation of the present invention is to construct Hankel matrices for each row and column of the two-dimensional magnetic resonance spectrum, and to obtain the Hankel matrix by using the non-convex h εThe norm fully constrains the low-rank characteristics of the Hankel matrix; the alternating direction multiplication method is then used to solve the established reconstruction model, and the optimal solution to the sub-problem of the Hankel matrix is ​​directly calculated based on this, which significantly improves the low-rank property and estimation accuracy of the obtained Hankel matrix. At the same time, the sub-problem of the magnetic resonance spectrum is quickly and approximately solved to obtain high-resolution reconstruction results.

[0035] The beneficial effects of the present invention are as follows: the Hankel matrix constructed by the rows and columns of the two-dimensional magnetic resonance spectrum is used as the processing object, and its binary index structure is fully utilized; a non-convex h ε The low-rank property of the norm-constrained Hankel matrix can impose more adequate low-rank constraints on the matrix than the existing low-rank regularization function; the reconstruction model is optimized based on the alternating direction multiplication method, and then the sub-problems about each variable are accurately solved, so that the obtained reconstructed magnetic resonance spectrum can restore a large amount of spectral peak information and remove artifacts, with higher resolution.

[0036] The present invention is mainly verified by simulation experiment method, and all steps and conclusions are verified to be correct on MATLAB8.0. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a flowchart of the workflow of the present invention;

[0038] Figure 2 is the fully sampled spatiotemporal correlation spectrum used in the simulation of the present invention;

[0039] Figure 3 It is the reconstruction result of the spatiotemporal correlation spectrum with a sampling rate of 25% using different methods;

[0040] Figure 3 It is the reconstruction result of the spatiotemporal correlation spectrum with a sampling rate of 25% using various methods (CS method, LRBHM method and the method of the present invention in turn);

[0041] Figure 4 It is the error of the reconstruction result of the spatiotemporal correlation spectrum with a sampling rate of 25% using various methods (CS method, LRBHM method and the method of the present invention in sequence). DETAILED DESCRIPTION

[0042] Reference Figure 1 The present invention is a magnetic resonance spectrum reconstruction method based on simultaneous low-rank constraints on rows and columns, and the specific steps include the following:

[0043] Step 1: construct a Hankel matrix for each row and column of the two-dimensional magnetic resonance spectrum.

[0044] (1a) Input a two-dimensional magnetic resonance spectroscopy signal's spatiotemporal encoding data y and an undersampling matrix U. For a two-dimensional magnetic resonance spectroscopy signal X to be reconstructed with N row vectors and M column vectors, the element in the mth column of the nth row is represented as x n，m ;

[0045] (1b) For the vector x in the mth column of X m , construct its corresponding Hankel matrix It can be expressed as:

[0046]

[0047] in, represents the linear operator that transforms a column vector into a Hankel matrix, and L represents the size parameter;

[0048] (1c) For the vector x in the nth row of X n , construct its corresponding Hankel matrix It can be expressed as:

[0049]

[0050] in, represents the linear operator that transforms a row vector into a Hankel matrix and K represents the size parameter.

[0051] Step 2: Based on the binary index structure of the two-dimensional magnetic resonance spectrum, the Hankel matrix constructed for each row and column is modified using h ε The norm constrains the low-rank characteristics to establish a magnetic resonance spectroscopy reconstruction model.

[0052] (2a) Using non-convex h ε The low-rank property of the Hankel matrix constructed with norm constraints, h represents any matrix L ε The norm is defined as:

[0053]

[0054] Among them, σ r is the rth singular value of L, ε>0 is a small constant;

[0055] (2b) Based on the binary index structure of two-dimensional magnetic resonance spectroscopy, h ε The norm simultaneously constrains the low-rank properties of the Hankel matrix constructed by rows and columns in X, thereby establishing a reconstruction model for two-dimensional magnetic resonance spectroscopy:

[0056]

[0057] Among them, F M represents the Fourier transform matrix of size M×M, IN represents the identity matrix of size N×N, represents the Kronecker product operator of matrices, vec(·) represents the operator that stacks the columns of a matrix into a vector, represents the square of the second norm of the vector, and λ represents the regularization parameter of the spatiotemporal encoding data fidelity term.

[0058] Step 3, solving the established reconstruction model by alternating direction multiplication method.

[0059] (3a) Introducing intermediate variables into the reconstruction model and The augmented Lagrangian function corresponding to the reconstruction model The expression is:

[0060]

[0061] Among them, β represents the penalty parameter, T m Indicates for L m The Lagrange multiplier introduced, Z n Indicates that for R n The introduced Lagrange multiplier, <·,·> represents the inner product operator of the matrix, The square of the Frobenius norm of the matrix is ​​expressed, and each optimization variable in the reconstruction model is solved according to the alternating direction multiplier method;

[0062] (3b) Fix X and R n , find the solution about L m Sub-questions:

[0063]

[0064] This subproblem can be solved by singular value threshold processing, which can be done by following the steps below:

[0065] (3b1) for the matrix Perform singular value decomposition and get:

[0066]

[0067] (3b2) Calculate the threshold parameter μ for the matrix singular values:

[0068]

[0069] (3b3) Perform threshold processing on all diagonal elements in the singular value matrix Δ:

[0070]

[0071] Where Th(·) represents the threshold function, and δ represents any diagonal element in Δ;

[0072] (3b4) The result of thresholding all diagonal elements in Δ can be expressed as Th(Δ), then m The optimal solution expression of the sub-problem is:

[0073] L m =PTh(Δ)V H

[0074] (3c) Fix X and L m , find the solution about R n Sub-questions:

[0075]

[0076] This sub-problem can also be solved by singular value thresholding, which can be done by following the steps below:

[0077] (3c1) for the matrix Perform singular value decomposition and get:

[0078]

[0079] Where E represents the left singular vector matrix, Q represents the right singular vector matrix, and Σ represents the singular value matrix;

[0080] (3c2) Threshold function Th(·) is used to perform threshold processing on all diagonal elements in the singular value matrix Σ. The result can be expressed as Th(Σ). Then, with respect to R n The optimal solution expression of the sub-problem is:

[0081] R n =ETh(Σ)Q H

[0082] (3d) Fixed L m and R n , solve the subproblem about X:

[0083]

[0084] This subproblem is a large-scale least squares problem, which can be quickly approximated by following the steps below:

[0085] (3d1) Definition and The conjugate operation of is defined as:

[0086]

[0087]

[0088] in, express The conjugate operator of , A represents any L×(Ν-L+1) size matrix, a k，i represents the element in the kth row and the ith column of A, express The conjugate operator of , B represents a matrix of any size K×(M-K+1), b k，i represents the element in the kth row and the ith column of B, min(·,·) represents the minimum function, max(·,·) represents the maximum function, (·) T represents the transpose operator of a matrix;

[0089] (3d2) Using the conjugate operator and Calculate the auxiliary matrices J1 and J2:

[0090]

[0091]

[0092] In addition, calculate the N×N size diagonal weight matrix With a diagonal weight matrix of size M×M in The i-th diagonal element of and The i-th diagonal element of The expression is:

[0093]

[0094]

[0095] (3d3) Introducing the intermediate variable s, the calculation formula for the intermediate variable s is:

[0096]

[0097] Among them, γ>0 is a small constant, I MN represents the identity matrix of size MN×MN, I M represents the identity matrix of size M×M, (·) -1 represents the inverse of the matrix, X last represents the solution of this subproblem in the previous iteration;

[0098] (3d4) After obtaining the intermediate variable s, the approximate solution expression of the subproblem about X is:

[0099]

[0100] in, is the Fourier transform matrix F of size M×MM The conjugate transposed matrix, U H It is the conjugate transposed matrix of the undersampling matrix U. Converting the obtained vec(X) into the matrix X is the approximate solution.

[0101] Step 4: After obtaining the solutions of the optimization variables in the reconstruction model, update the Lagrange multiplier T m and Z n , and the penalty parameter β:

[0102]

[0103]

[0104] β=τβ

[0105] Where τ>1 represents a growth factor, and steps (3) to (4) are repeated until the reconstructed magnetic resonance spectrum meets the conditions or the number of iterations reaches a preset upper limit.

[0106] The effect of the present invention can be further illustrated by the following simulation experiment:

[0107] 1. Experimental conditions and contents

[0108] Experimental conditions: The experiment used a Poisson gap sampling template with a sampling rate of 20%; the experimental magnetic resonance spectrum used a real time-space correlation spectrum such as Figure 2 As shown in the figure, the signal-to-noise ratio (SNR) is used as an objective evaluation index of the reconstruction result, and its definition is:

[0109]

[0110] Where X is the real magnetic resonance spectrum, To reconstruct the magnetic resonance spectrum. The higher the SNR value, the closer the reconstructed magnetic resonance spectrum is to the real magnetic resonance spectrum.

[0111] Experimental content: Under the above experimental conditions, the representative CS method and LRBHM method in the field of magnetic resonance spectroscopy reconstruction were selected for comparison with the method of the present invention.

[0112] Experiment 1: The method of the present invention, the CS method and the LRBHM method were used to analyze Figure 2 The CS method is based on the sparsity of constrained magnetic resonance spectra and uses an iterative soft threshold method to solve the reconstruction model. The reconstruction result is Figure 3 The reconstruction error is Figure 4 The LRBHM method first constructs the Hankel matrix for each row of the two-dimensional magnetic resonance spectrum, and then uses these Hankel matrices to construct a larger block Hankel matrix, and uses the nuclear norm to constrain its low-rank characteristics. The reconstruction result is Figure 3 The reconstruction error is Figure 4 b; In the experiment, the method of the present invention sets the size parameter L to be an integer of N / 2, the size parameter K to be an integer of M / 2, and h ε The small constant ε in the norm is 10 -3 , the regularization parameter of the spatiotemporal coding data fidelity term is λ=10 5 , the approximate solution of the subproblem of X has a small constant γ = 0.01, a growth factor τ = 1.05, and the reconstruction result of the method of the present invention is Figure 3 The reconstruction error is c. Figure 4 c.

[0113] from Figure 3 It can be seen from the reconstruction results that the CS method causes serious artifacts in the reconstructed magnetic resonance spectrum and only recovers a small number of weakened spectral peaks; although the LRBHM method recovers all spectral peaks, it causes the range of some spectral peaks to shrink and loses a certain amount of spectral peak intensity; the method of the present invention not only effectively suppresses artifacts, but also retains most of the spectral peak information in all spectral peaks. Figure 4 It can be seen from the reconstruction error diagram that the reconstruction error of the method of the present invention is significantly less than that of other methods, so the corresponding reconstruction result is closest to the real magnetic resonance spectrum.

[0114] Table 1 SNR indicators of different magnetic resonance spectroscopy reconstruction methods

[0115] Magnetic resonance spectroscopy CS Method LRBHM Method Method of the present invention Space-time correlation spectroscopy 8.62 22.03 34.72

[0116] Table 1 gives Figure 3 From the SNR indicators of the reconstruction results of each method, it can be seen that the SNR value of the method of the present invention is significantly improved compared with other methods, indicating that the quality of the reconstructed magnetic resonance spectrum obtained by the method of the present invention is the best. This result is consistent with the reconstruction effect diagram.

[0117] The above experiments show that the reconstructed magnetic resonance spectrum obtained by the method of the present invention visually eliminates a large number of artifacts and restores sufficient spectral peak information, and the corresponding objective evaluation index is significantly higher than that of the comparison method. This shows that the present invention is effective for magnetic resonance spectrum reconstruction.

Claims

1. A method for magnetic resonance spectroscopy reconstruction based on simultaneous low-rank constraints of rows and columns, comprising the following steps: (1) Input a two-dimensional magnetic resonance spectroscopy signal's spatiotemporal encoding data y and an undersampling matrix U. For a two-dimensional magnetic resonance spectroscopy signal X to be reconstructed with N row vectors and M column vectors, the element in the mth column of the nth row is represented as x n，m ; (2) For the vector x in the mth column of X m , construct its corresponding Hankel matrix It can be expressed as: in, represents the linear operator that transforms a column vector into a Hankel matrix, and L represents the size parameter; (3) For the vector x in the nth row of X n , construct its corresponding Hankel matrix It can be expressed as: in, represents the linear operator that transforms a row vector into a Hankel matrix, and K represents the size parameter; (4) Using non-convex h ε The low-rank property of the Hankel matrix constructed with norm constraints, h represents any matrix L ε The norm is defined as: Among them, σ r is the rth singular value of L, ε>0 is a small constant; (5) When using h ε Based on the low-rank characteristics of the Hankel matrix constructed by simultaneously constraining the rows and columns in X, a reconstruction model for two-dimensional magnetic resonance spectroscopy is established: Among them, F M represents the Fourier transform matrix of size M×M, I N represents the identity matrix of size N×N, represents the Kronecker product operator of matrices, vec(·) represents the operator that stacks the columns of a matrix into a vector, represents the square of the second norm of the vector, λ represents the regularization parameter of the spatiotemporal coding data fidelity term; (6) Introducing intermediate variables into the reconstruction model and And the sub-problems of each optimization variable in the reconstruction model are solved according to the alternating direction multiplier method: (6a) Fix X and R n , find the solution about L m Sub-questions: Among them, β represents the penalty parameter, T m Indicates for L m The Lagrange multiplier introduced, Represents the square of the Frobenius norm of the matrix; (6b) Fix X and L m , find the solution about R n Sub-questions: Among them, Z n Indicates that for R n The Lagrange multiplier introduced; (6c) Fixed L m and R n , solve the subproblem about X: (7) After obtaining the solution of each optimization variable, update the Lagrange multiplier T m and Z n , and the penalty parameter β: β=τβ Where τ>1 represents a growth factor, and steps (6) to (7) are repeated until the reconstructed magnetic resonance spectrum meets the conditions or the number of iterations reaches a preset upper limit.

2. The magnetic resonance spectroscopy reconstruction method based on simultaneous row and column low-rank constraints according to claim 1, characterized in that: Steps (6a) and (6b) solve for L m and R n The sub-problem can be solved by following the steps below: (6a1) For the matrix Perform singular value decomposition and get: Where P represents the left singular vector matrix, V represents the right singular vector matrix, and Δ represents the singular value matrix; (6a2) Calculate the threshold parameter μ for the matrix singular values: (6a3) Perform threshold processing on all diagonal elements in the singular value matrix Δ: Where Th(·) represents the threshold function, and δ represents any diagonal element in Δ; (6a4) The result of thresholding all diagonal elements in Δ can be expressed as Th(Δ), then m The optimal solution expression of the sub-problem is: L m =PTh(Δ)V H (6b1) for the matrix Perform singular value decomposition and get: Where E represents the left singular vector matrix, Q represents the right singular vector matrix, and Σ represents the singular value matrix; (6b2) Threshold function Th(·) is used to perform threshold processing on all diagonal elements in the singular value matrix Σ. The result can be expressed as Th(Σ). Then, with respect to R n The optimal solution expression of the sub-problem is: R n =ETh(Σ) Q H 。 3. The magnetic resonance spectrum reconstruction method based on simultaneous row and column low-rank constraints according to claim 1, characterized in that: To solve the subproblem about X in step (6c), we can follow these steps: (6c1) This subproblem is a large-scale least squares problem. To quickly calculate its approximate solution, we need to define and The conjugate operation of is defined as: in, express The conjugate operator of , A represents any L×(Ν-L+1) size matrix, a k，i represents the element in the kth row and the ith column of A, express The conjugate operator of , B represents a matrix of any size K×(M-K+1), b k，i represents the element in the kth row and the ith column of B, min(·,·) represents the minimum function, max(·,·) represents the maximum function, (·) T represents the transpose operator of a matrix; (6c2) Using the conjugate operator and Calculate the auxiliary matrices J1 and J2: In addition, calculate the N×N size diagonal weight matrix With a diagonal weight matrix of size M×M in The i-th diagonal element of and The i-th diagonal element of The expression is: (6c3) Introduce the intermediate variable s, and the calculation formula for the intermediate variable s is: Among them, γ>0 is a small constant, I MN represents the identity matrix of size MN×MN, I M represents the identity matrix of size M×M, I N represents the identity matrix of size N×N, represents the Kronecker product operator of matrices, (·) -1 represents the inverse of the matrix, X last represents the solution of this subproblem in the previous iteration, and vec(·) represents the operator that stacks the columns of the matrix into a vector; (6c4) After obtaining the intermediate variable s, the approximate solution expression of the subproblem about X is: in, is the Fourier transform matrix F of size M×M M The conjugate transposed matrix, U H It is the conjugate transposed matrix of the undersampling matrix U. Converting the obtained vec(X) into the matrix X is the approximate solution.

Images