A magnetic resonance spectrum denoising method based on exponential decomposition constraints

Abstract

A magnetic resonance spectrum denoising method based on exponential decomposition constraints, relating to a magnetic resonance spectrum denoising method. The method comprises the following steps: 1) proposing a magnetic resonance spectrum denoising model that simultaneously utilizes the low-rank characteristics of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint; 2) denoising the magnetic resonance spectrum through an iterative algorithm. The method is a magnetic resonance spectrum denoising method that simultaneously utilizes the low-rank characteristics of the Hankel matrix of the magnetic resonance time domain signal and the characteristics of being decomposed into a plurality of exponential signal superposition characteristics, and has the characteristics of high noise reduction quality.

Description

Technical Field

[0001] The invention relates to a magnetic resonance spectrum denoising method, in particular to a method for solving the magnetic resonance spectrum denoising problem based on exponential decomposition constraints. Background Art

[0002] As we all know, magnetic resonance spectroscopy is one of the important analytical tools in the fields of biology, chemistry, and medicine. However, the acquisition process of magnetic resonance spectroscopy is often affected by noise, and there is a problem of low signal-to-noise ratio. In practical applications, repeated acquisition and averaging are often used to reduce noise. However, multiple repeated acquisitions will increase the sampling time and cause a waste of time. Magnetic resonance spectroscopy denoising by signal processing is to extract a noise-free signal from the original noisy signal, usually modeled with Gaussian noise, thereby improving the signal-to-noise ratio of the nuclear magnetic resonance signal.

[0003] In the past few decades, researchers have proposed many denoising methods for magnetic resonance spectra. For example, Cadzow (Yung, Ya Lin, Lian Pin, Hwang, "NMR signal enhancement based on matrix property mappings," Journal of Magnetic Resonance, Series A, 103, 109-114, 1993) is a classic denoising method, which denoises by specifying the Hankel matrix structure of the truncated matrix, but the solution obtained is not necessarily the global optimal solution. In magnetic resonance spectrum denoising, many researchers use the mathematical properties of magnetic resonance signals to denoise and reconstruct magnetic resonance spectra, such as the low rank properties of the Hankel matrix of the magnetic resonance spectrum time domain signal (Xiaobo Qu, Maxim, Mayzel, Jian-Feng Cai, Zhong Chen, Vladislav Orekhov, "Accelerated NMR spectroscopy with low-rank reconstruction," Angewandte Chemie International Edition, 54, 852-854, 2015).The method of Qu et al. denoises the signal by constructing a Lagrangian convex optimization solution model based on the Hankel matrix nuclear norm, and applies it to higher dimensional spectral reconstruction (Jiaxi Ying, Hengfa Lu, Qingtao Wei, Jian-Feng Cai, Di Guo, Jihui Wu, Zhong Chen, Xiaobo Qu, "Hankel matrix nuclearnorm regularized tensor completion for N-dimensional exponential signals," IEEE Transactions on Signal Processing, 65(14), 3702-3717, 2017) and time-domain and frequency-domain mixed signal reconstruction (Hengfa Lu, Xinlin Zhang, Tianyu Qiu, Jian Yang, Jiaxi Ying, Di Guo, Zhong Chen, Xiaobo Qu, "Low rank enhanced matrix recovery of hybrid time andfrequency data in fast magnetic resonance spectroscopy," IEEE Transactions on Biomedical Engineering, 65(4), 809-820, 2018). Recently, Ying et al. proposed a method combining Vandermonde decomposition to optimize the original algorithm (Jiaxi Ying, Jian-Feng Cai, Di Guo, Gongguo Tang, ZhongChen, Xiaobo Qu, "Vandermonde factorization of Hankel matrix for complex exponential signal recovery-application infast NMR spectroscopy," IEEE Transactions on Signal Processing, 66(21), 5520-5533, 2018), making the low-rank constraint more reasonable. However, this method takes a long time to denoise when the number of sampling points is large.(Tianyu Qiu, Wenjing Liao, Yihui Huang, Jinyu Wu, Di Guo, Dongbao Liu, Xinwang, Jian-Feng Cai, Bingwen Hu, Xiaobo Qu, "An automatic denoising method for NMR spectroscopy based on low-rank Hankel model," IEEE Transactions on Instrumentation and Measurement, 70, 1-12, 2021) A denoising method with automatic parameter selection is proposed to autonomously obtain the optimal solution. (Xiaobo Qu, Yihui Huang, Hengfa Lu, Tianyu Qiu, Di Guo, Vladislav Orekhov, Zhong Chen, "Accelerated Nuclear Magnetic Resonance Spectroscopy with Deep Learning," Angewandte Chemie International Edition, 59(26), 10297-10300, 2020) uses an exponential model to construct simulated signals to train neural networks, which can quickly obtain high-fidelity magnetic resonance spectrum reconstruction results. (Yihui Huang#, Jinkui Zhao#, ZiWang, Di Guo, Xiaobo Qu, "Exponential signal recovery with deep Hankel matrix factorization," IEEE Transactions on Neural Networks and Learning Systems, 1-13, 2021) A deep learning method using low-rank Hankel matrix factorization is proposed, and its neural network structure is constructed by imitating the iterative process of the exponential signal reconstruction method. This method can obtain lower reconstruction error and can effectively retain weak peaks.

[0004] In magnetic resonance spectroscopy denoising, the denoising quality of existing methods still needs to be improved, and the properties of the singular values ​​of the Hankel matrix have not been fully utilized. Summary of the invention

[0005] The object of the present invention is to provide a magnetic resonance spectrum denoising method based on exponential decomposition constraints, which can simultaneously utilize the low-rank characteristics of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint to achieve high-quality magnetic resonance spectrum denoising.

[0006] The present invention comprises the following steps:

[0007] 1) A magnetic resonance spectrum denoising model is proposed that simultaneously utilizes the low-rank property of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint;

[0008] 2) De-noising the magnetic resonance spectrum through an iterative algorithm.

[0009] In step 1), the specific method of proposing a magnetic resonance spectrum denoising model that simultaneously utilizes the low-rank characteristics of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint is:

[0010] use It represents the noise-free magnetic resonance time domain signal and can be expressed as Where N represents the length of the signal; Modeled as a superposition of decaying exponential signals, where The nth element x n It can be expressed as:

[0011]

[0012] in, n=0,...,N-1.c r 、f r and τ r are the amplitude, resonance frequency and attenuation coefficient of the rth peak respectively, r=1, 2, ..., R, where R represents the total number of signal spectrum peaks, and Δt represents the sampling interval.

[0013] Defining Operators P = N - Q + 1, Indicates that the FID signal The process of converting to the Hankel matrix, the Hankel matrix can be expressed as:

[0014]

[0015] When N is an odd number, Q=(N+1) / 2; when N is an even number, Q=N / 2.

[0016] The rank of the Hankel matrix of the FID signal is the same as the number of exponential signals (that is, the number of spectral peaks); when the number of exponentials R is less than the rank of the Hankel matrix When the dimension It has low-rank characteristics;

[0017] Suppose the index set of all peaks is |·| represents the number of elements in the index set. The magnitude of the index in formula (1) can be further expressed as And c r Arrange as vector Signal All z in r Can be arranged into a Vandermonde matrix in Will Expressed as

[0018] set up For the signal to be denoised, the Hankel matrix converted from x satisfies the low-rank property and has a decomposition form of x=Zc, that is, it has an exponential decomposition constraint.

[0019] The proposed denoising model is:

[0020]

[0021] in, is the measurement signal with noise, is the signal to be denoised, is a Vandermonde matrix. ||·|| * is the nuclear norm, which is a convex relaxation of the low-rank constraint, ||·||2 represents the 2-norm of the vector, and η represents the range of data fidelity.

[0022] This model can be further expressed as:

[0023]

[0024] Among them, λ and μ are regularization parameters, both of which are positive real numbers.

[0025] In step 2), the specific method of denoising the magnetic resonance spectrum by an iterative algorithm may be:

[0026] Formula (4) can be calculated by the alternating multiplier method; first, the augmented Lagrangian form of formula (4) is:

[0027]

[0028] Where D is the dual variable, <·,·> represents the inner product, represents the square of the norm of the matrix F, λ>0, β>0 and μ>0 are all regularization parameters.

[0029] The sub-problem of model (5) can be expressed as:

[0030]

[0031] The sub-problems of formula (6) are solved by the following method:

[0032] (1)xi+1 It can be obtained by solving the following formula:

[0033]

[0034] Taking the derivative of the above formula, we can get the solution as

[0035]

[0036] in, for The adjoint operator of .

[0037] (2) Z i+1 It can be obtained by solving the following formula:

[0038]

[0039] Therefore, we can get Z i+1 The steps to solve are:

[0040] Step 1: Calculate the Hankel matrix Singular value decomposition of in(·) H Represents the conjugate transpose of the matrix and then decomposes the matrix U∑V H Truncate to The subscript R means taking the first R columns of the original matrix;

[0041] Step 2: Calculate U R l M R ≈U R f The least squares solution M R , the superscripts l and f represent the deletion matrix U R The first row and deletion of the matrix U R The last line of . Calculate M R , and regard it as

[0042] Step 3: r=1,2,...,R is arranged into a Vandermonde matrix Z i+1 ,

[0043] (1)c i+1 It can be obtained by taking the derivative of the following formula

[0044]

[0045] Its solution is

[0046] c i+1 =((Z i+1 ) H Zi+1 ) -1 (Z i+1 ) H x i+1 (11)

[0047] in,(·) -1 It means to find the pseudo-inverse.

[0048] (2)U i+1 It can be obtained by solving the following formula:

[0049]

[0050] The solution is:

[0051]

[0052] (3)D i+1 The solution is:

[0053]

[0054] Where t represents the update step size.

[0055] The entire iterative process is:

[0056] Initialization: input y, As well as the total number of spectral peaks R, set β = t = 1, μ = 0.1, and the range of data fidelity η = 10 -8 . Dual variable D 0 =1, variable x 0 =y,Z 0 =0, c 0 =0,U 0 =0, Δx=10 -10 , where 1 and 0 represent matrices or vectors with all 1s and all 0s, respectively, of the same dimension as the corresponding matrix.

[0057] When Δx>η, repeat the calculation steps a)-e) and calculate the error of the FID of two adjacent de-noised magnetic resonance spectra. When the FID error of two adjacent de-noised magnetic resonance spectra is When it is reduced to the specified range, that is, when Δx < η, the iterative operation is stopped, and the denoising result x of the last iteration is taken as the final denoising result.

[0058] The method designs a denoising method for magnetic resonance spectra, in particular, simultaneously utilizing the translation invariance of the singular value decomposition of the Hankel matrix and its low-rank characteristics. The method comprises: 1) proposing a magnetic resonance spectrum denoising model that simultaneously utilizes the low-rank characteristics of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint; 2) denoising the magnetic resonance spectrum by iterative operation. The present invention is used for denoising magnetic resonance spectra and has the characteristics of high denoising quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a noise-free spectrum in the embodiment of the present invention.

[0060] Figure 2 is the noisy spectrum input in the embodiment of the present invention.

[0061] Figure 3 It is the spectrum after the noisy spectrum in the embodiment of the present invention is denoised. DETAILED DESCRIPTION

[0062] To make the technical solution of the present invention clearer, the following embodiments will further illustrate the present invention in conjunction with the accompanying drawings. It should be understood that the preferred embodiments are only for illustrating the present invention, rather than for limiting the protection scope of the present invention.

[0063] The embodiment of the present invention uses an exponential function to generate a noise-free magnetic resonance time domain signal, the signal length is N = 512, and the total number of spectrum peaks is R = 5. By performing Fourier transform on it, a noise-free spectrum can be obtained, such as Figure 1 As shown, a Gaussian white noise with a standard deviation of 0.02 is added to the time domain signal, and a Fourier transform is performed to obtain the spectrum with the noise added, as shown in Figure 2 shown.

[0064] The embodiment of the present invention includes the following steps:

[0065] 1) A magnetic resonance spectrum denoising model is proposed that simultaneously utilizes the low-rank property of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint:

[0066] use represents the noise-free magnetic resonance time domain signal, expressed as Where N represents the length of the signal; Modeled as a superposition of decaying exponential signals, where The nth element x n It is expressed as:

[0067]

[0068] in, n=0,...,N-1.c r 、f r and τr are the amplitude, resonance frequency and attenuation coefficient of the rth peak respectively, r=1, 2, ..., R, where R represents the total number of signal spectrum peaks, and Δt represents the sampling interval.

[0069] Defining Operators P = N - Q + 1, Indicates that the FID signal The process of converting to the Hankel matrix, the Hankel matrix can be expressed as

[0070]

[0071] When N is an odd number, Q=(N+1) / 2; when N is an even number, Q=N / 2.

[0072] The rank of the Hankel matrix of the FID signal is the same as the number of exponential signals (i.e., the number of spectral peaks). When the dimension It has low-rank properties.

[0073] Suppose the index set of all peaks is |·| represents the number of elements in the index set. The magnitude of the index in formula (1) can be further expressed as And c r Arrange as vector Signal All z in r Can be arranged into a Vandermonde matrix in Will Expressed as

[0074] set up For the signal to be denoised, the Hankel matrix converted from x satisfies the low-rank property and has a decomposition form of x=Zc, that is, it has an exponential decomposition constraint.

[0075] The proposed denoising model is:

[0076]

[0077] in, is the measurement signal with noise, is the signal to be denoised, is a Vandermonde matrix. ||·|| * is the nuclear norm, which is a convex relaxation of the low-rank constraint, ||·||2 represents the 2-norm of the vector, and η represents the range of data fidelity.

[0078] This model can be further expressed as:

[0079]

[0080] Among them, λ and μ are regularization parameters, both of which are positive real numbers.

[0081] 2) De-noising the magnetic resonance spectrum by iterative algorithm:

[0082] Model (4) can be calculated by the alternating multiplier method. First, the augmented Lagrangian form of (4) is:

[0083]

[0084] Where D is the dual variable, <·,·> represents the inner product, represents the square of the norm of the matrix F, λ>0, β>0 and μ>0 are all regularization parameters.

[0085] The sub-problem of model (5) can be expressed as:

[0086]

[0087] The sub-problems of formula (6) are solved by the following method:

[0088] a)x i+1 It can be obtained by solving the following formula

[0089]

[0090] Taking the derivative of the above formula, we can get the solution as

[0091]

[0092] in, for The adjoint operator of .

[0093] b) Z i+1 It can be obtained by solving the following formula

[0094]

[0095] Therefore, we can get Z i+1 The solution steps are:

[0096] Step 1: Calculate the Hankel matrix Singular value decomposition of in(·) H Represents the conjugate transpose of the matrix and then decomposes the matrix U∑V H Truncate to The subscript R means taking the first R columns of the original matrix;

[0097] Step 2: Calculate U R l M R ≈U R f The least squares solution M R , the superscripts l and f represent the deletion matrix U R The first row and deletion of the matrix U R The last line of . Calculate M R , and regard it as

[0098] Step 3: Arranged into the Vandermonde matrix Z i+1 ,

[0099] c)c i+1 It can be obtained by taking the derivative of the following formula:

[0100]

[0101] The solution is:

[0102] c i+1 =((Z i+1 ) H Z i+1 ) -1 (Z i+1 ) H x i+1 (11)

[0103] in,(·) -1 It means to find the pseudo-inverse.

[0104] d)U i+1 It can be obtained by solving the following formula:

[0105]

[0106] The solution is:

[0107]

[0108] e)D i+1 The solution is

[0109]

[0110] Where t represents the update step size.

[0111] The entire iterative process is:

[0112] Initialization: input y, As well as the total number of spectral peaks R, set β = t = 1, μ = 0.1, and the range of data fidelity η = 10 -8 . Dual variable D 0 =1, variable x 0 =y,Z 0 =0, c 0 =0,U 0 =0, Δx=10 -10 , where 1 and 0 represent matrices or vectors with all 1s and all 0s, respectively, of the same dimension as the corresponding matrix.

[0113] When Δx>η, repeat the calculation steps a)-e) and calculate the error of the FID of two adjacent de-noised magnetic resonance spectra. When the FID error of two adjacent de-noised magnetic resonance spectra is When it is reduced to the specified range, that is, when Δx < η, the iterative operation is stopped, and the denoising result x of the last iteration is taken as the final denoising result.

[0114] For the denoised time domain signal x out By performing Fourier transform, we can get the denoised spectrum, such as Figure 3 shown.

Claims

1. A magnetic resonance spectrum denoising method based on exponential decomposition constraints, It is characterized by the following steps: 1) A magnetic resonance spectrum denoising model is proposed that simultaneously utilizes the low-rank characteristics of the Hankel matrix of the magnetic resonance spectrum time domain signal and the exponential decomposition constraint. The specific method is as follows: use represents the noise-free magnetic resonance time domain signal, expressed as in, N represents the length of the signal; Modeled as a superposition of decaying exponential signals, where The nth element x n It is expressed as: in, n=0,...,N-1;c r 、f r and τ r are the amplitude, resonance frequency and attenuation coefficient of the rth peak respectively; r = 1, 2, ..., R, where R represents the total number of signal spectrum peaks and Δt represents the sampling interval; Defining Operators P=N-Q+1, R represents the FID signal The process of converting to the Hankel matrix, the Hankel matrix is ​​expressed as: When N is an odd number, Q = (N + 1) / 2; when N is an even number, Q = N / 2; The rank of the Hankel matrix of the FID signal is the same as the number of exponential signals; when the number of exponential signals R is less than the rank of the Hankel matrix When the dimension It has low-rank characteristics; Suppose the index set of all peaks is |·| represents the number of elements in the index set, and the magnitude of the exponent in formula (1) is expressed as And c r Arrange as vector Signal All z in r Arranged into a Vandermonde matrix in Will Expressed as set up is the signal to be denoised, the Hankel matrix converted from x satisfies the low-rank property and has the decomposition form of x = Zc, that is, it has an exponential decomposition constraint; The proposed denoising model is: in, is the measurement signal with noise, is the signal to be denoised, is a Vandermonde matrix; ||·|| * is the nuclear norm, which is a convex relaxation of the low-rank constraint, ||·||2 represents the 2-norm of the vector, and η represents the range of data fidelity; The model is represented as: Among them, λ and μ are regularization parameters, both of which are positive real numbers; 2) De-noising the magnetic resonance spectrum by an iterative algorithm. The specific method is as follows: Formula (4) is calculated by the alternating multiplier method; first, the augmented Lagrangian form of formula (4) is: Where D is the dual variable, <·,·> represents the inner product, represents the square of the norm of the matrix F, λ>0, β>0 and μ>0 are all regularization parameters; The sub-problem of formula (5) is expressed as: The sub-problems of formula (6) are solved by the following method: (1)x i+1 By solving the following formula: Taking the derivative of the above formula, the solution is: in, for The adjoint operator of ; (2) Z i+1 By solving the following formula: So we get Z i+1 The steps to solve are: Step 1: Calculate the Hankel matrix Singular value decomposition of in(·) H Represents the conjugate transpose of the matrix and then decomposes the matrix U∑V H Truncate to The subscript R means taking the first R columns of the original matrix; Step 2: Calculate U R l M R ≈U R f The least squares solution M R , the superscripts l and f represent the deletion matrix U R The first row and deletion of the matrix U R The last line of ; calculate M R , and regard it as Step 3: Arranged into a Vandermonde matrix (3)c i+1 By taking the derivative of the following formula: The solution is: c i+1 =((Z i+1 ) H z i+1 ) -1 (Z i+1 ) H x i+1 (11) in,(·) -1 It means to find the pseudo-inverse; (1)U i+1 By solving the following formula: The solution is: (2)D i+1 The solution is: Where t represents the update step size; The entire iterative process is: Initialization: input y, As well as the total number of spectral peaks R, set β = t = 1, μ = 0.1, and the range of data fidelity η = 10 -8 ; Dual variable D 0 =1, variable x 0 =y,Z 0 =0, c 0 =0,U 0 =0, Δx=10 -10 , where 1 and 0 represent matrices or vectors with all 1s and all 0s, respectively, of the same dimension as the corresponding matrix; When Δx>η, repeat steps a)-e) and calculate the error of the FID of two adjacent de-noised magnetic resonance spectra. When the FID error of two adjacent de-noised magnetic resonance spectra is When it is reduced to the specified range, that is, Δx<η, the iterative operation is stopped, and the denoising result x of the last iteration is used as the final denoising result.

Images