An improved ESPIRiT reconstruction method based on sparse transform learning

Abstract

The present invention relates to an improved ESPIRiT reconstruction method based on sparse transform learning, and belongs to the technical field of magnetic resonance imaging. The present invention introduces data-driven adaptive sparse transform learning (Transform Learning, TL) into the ESPIRiT model, and proposes a parallel MRI reconstruction algorithm based on ESPIRiT multiple groups of sensitivity operators combined with data-driven adaptive sparse transform learning regularization terms, named TL-ESPIRiT. The algorithm introduces tight frame technology for solution, and simulation experimental results show that compared with the LpJTV-ESPIRiT and TV-ESPIRiT models, the TL-ESPIRiT model reduces the staircase artifacts of the reconstructed image, the texture details and edge contour information are kept intact, and the reconstructed image is clearer.

Description

Technical Field

[0001] The invention relates to an improved ESPIRiT reconstruction method based on sparse transformation learning, and belongs to the technical field of magnetic resonance imaging. Background Art

[0002] Magnetic resonance imaging (MRI) provides a variety of contrast mechanisms, good visualization of physiological functions, and has the advantages of being non-invasive, non-destructive, and rich in image information, making it one of the main imaging diagnostic methods. However, MRI data acquisition takes a long time, which leads to longer examination times and higher imaging costs. Therefore, the main problem that magnetic resonance imaging technology solves is to speed up magnetic resonance imaging scanning and improve reconstruction quality.

[0003] Parallel imaging (PI) is an imaging technology based on multi-channel coils. It uses a nuclear magnetic resonance phased array to obtain information from multiple closely positioned nuclear magnetic resonance receiving coils, reducing the number of phase encoding steps required for imaging and speeding up scanning. In 2000, Griswold et al. proposed parallel imaging of local sensitivity (PILS). The spatial encoding of this method is performed in the RF coil array to speed up image acquisition. In 2002, Griswold et al. combined the advantages of PILS and VD-AUTO-SNASH and proposed the generalized automatic calibration partial parallel acquisition (GRAPPA) algorithm. Although the introduction of parallel imaging has further developed the rapid imaging technology of magnetic resonance, the scanning speed is limited and the artifacts of the reconstructed image cannot be eliminated well.

[0004] Compressed sensing theory shows that if the signal is sparse, the original signal can be accurately restored from a small amount of measured data or undersampled data through an optimization algorithm. This process is called sparse recovery. In 2009, Dong et al. proposed a CS-SENSE algorithm, which combines sensitivity encoding (SENSE) and fast magnetic resonance imaging compressed sensing (SparseMRI) to improve the compression factor. In the same year, Lustig et al. proposed an automatic calibration parallel imaging technology based on compressed sensing, which is used for the corrected sparse model of multi-coil images. In 2014, Uecker et al. combined the advantages of SENSE and GRAPPA and proposed an iterative self-consistent parallel imaging reconstruction using eigenvector operators (ESPIRiT: iTerative Self-consistentParallel Imaging Reconstruction using Eigenvector maps) algorithm. The algorithm achieves coil image reconstruction by using the calibration kernel filter of K space and the pre-known sensitivity operator information. In 2019, Duan et al. introduced the Lp-norm joint total variation regularization term (LpJTV) into the ESPIRiT model, transforming the reconstruction problem into a minimization problem with composite terms. Using the Bregman iteration technique, a series of unconstrained problems were obtained, and the operator splitting technique was used to transform the unconstrained problem into a gradient problem and a denoising problem for solution. The introduction of the regularization term enables the ESPIRiT model to show good performance in suppressing noise and removing blur artifacts.

[0005] Parallel MRI reconstruction based on the ESPIRiT model requires regularization terms to suppress noise and remove blurring artifacts. The ESPIRiT model is a parallel imaging reconstruction method that uses iterative self-consistency of eigenvector operators, which has strong flexibility and compatibility in MRI reconstruction. Currently, the common regularizer-based ESPIRiT models mainly include the TV-ESPIRiT algorithm based on the TV regularization term and the LpJTV-ESPIRiT algorithm based on the Lp pseudo-norm joint total variation regularization term. Although these algorithms can reconstruct images, the reconstructed images have staircase artifacts and blurring artifacts, and the details and edges are less intact. Summary of the invention

[0006] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a parallel magnetic resonance imaging reconstruction method based on sparse transform learning and ESPIRiT, which can further improve the quality of magnetic resonance imaging reconstructed images.

[0007] The object of the present invention is achieved through the following technical solution: a parallel magnetic resonance imaging reconstruction method based on sparse transform learning and ESPIRiT, which comprises the following steps:

[0008] S0: Initialize, set k = 0, x 0 =(RFS) H y, t 0 = 0, b 0 = H(Φ 0 z 0 , λ),

[0009] where k is the loop variable, representing the k-th iteration of the data, and the superscript " 0 " represents the initial value. represents the multi-component image to be reconstructed, N = N h ×N v represents the number of pixels of the single-coil magnitude image, N v and N h are the number of rows and columns of the single-coil magnitude image respectively, x j represents the j-th component image to be reconstructed, j = 1......J represents the index of the sensitivity map group, x 1 and x J represent the first component image and the J-th component image of the column-vectorized multi-component image x to be reconstructed respectively, J represents the total number of groups of sensitivity maps, and J takes 2 in the present invention, x 0 represents the initial value of x. represents the undersampled multi-coil k-space data, represents the undersampled k-space data of the c-th coil, M is the actual number of sampled points of the single-coil k-space data, and M << N, c = 1......C represents the index of the receiving coils, y 1 and y C represent the first coil data and the C-th coil data of the column-vectorized undersampled multi-coil k-space data y respectively, C represents the number of receiving coils used for parallel imaging, and the superscript " T " represents matrix transpose, and the superscript " H " represents matrix complex conjugate transpose. is an undersampling matrix, where represents the matrix for selecting the sampling point positions from the k-space grid, I C is the C×C identity matrix, represents the Kronecker product. is the Fourier operator, and represent the two-dimensional Fourier transform matrices of N h points and N v points respectively. is the sensitivity map group matrix, obtained by the ESPIRiT model, and the j-th group sensitivity map of the c-th coil of S is expressed by S cj Indicates that S 11 The first group of sensitivity diagrams for the first coil, S CJ Represents the Jth group of sensitivity diagrams for the Cth coil. represents the intermediate variable, t 0 and are t and The corresponding initial value.

[0010] It represents the linear operator extracted from the i-th image block, which can extract J images of size The extracted J image blocks are converted into The matrix, P 1 (·) indicates the first matrix extraction linear operator, Indicates the Nth P The linear operator is extracted from the matrix. express The sparse transformation matrix of the image patch, W 0 represents the initial value of W, express The discrete cosine transform matrix of the point. Ψ is the defined transformation, and Ψ 0 is the initial value of Ψ; Φ is a tight frame, and Φ H Φ=I 2nN , Φ 0 is the initial value of Φ; H(x, θ) is the hard threshold function, x represents the input matrix, θ represents the threshold, λ is the parameter, and λ>0. z is the intermediate variable, and all z are concatenated into a matrix in order. And Z represents 2N image blocks extracted from z and spliced ​​into an n×2N matrix, P 1 z 1 and P N z 2 They represent the first image block and the 2Nth image block extracted from the multi-component image z by the intermediate variables, respectively. It means Z=(P 1 z 1 … P N z 1 P 1 z 2 … P N z 2 )'s initial value. Auxiliary variables is the column-vectorized encoding matrix, b 0 Indicates the initial value of b.

[0011] S1: Calculate the intermediate variable z for the k+1th iteration k+1 , the calculation formula is as follows:

[0012]

[0013] in, Represents the intermediate variable of the kth iteration L is The Lipschitz constant of the gradient.

[0014] S2: Initialize j=1;

[0015] S3: Update B k = unvec(b k ), the calculation formula is as follows:

[0016]

[0017] S4: Calculate the k+1th intermediate variable, Z k+1 :

[0018]

[0019] S5: Z k+1 (B k ) H Perform singular value decomposition and get Z k+1 (Bk) H =U∑V H ;

[0020] Among them, B k is the kth iteration of the auxiliary variable B, Z k+1 is the k+1th iteration of the intermediate variable Z, and U, ∑, V are the singular value decomposition factors.

[0021] S6: Transform update, calculate the sparse transformation W of the k+1th iteration k+1 , the calculation formula is as follows:

[0022] W k+ 1=VU H .

[0023] S7: Hard threshold denoising, calculate the auxiliary variable b for the k+1th iteration k+1 , the calculation formula is as follows:

[0024] b k+1 =H(Φ k z k ,λ).

[0025] S8: Determine whether all image components are completed. If j=J, go to step S9; otherwise, set j=j+1 and return to S3.

[0026] S9: Calculate the image x to be reconstructed in the (k + 1)-th iteration k+1 , and the calculation formula is as follows:

[0027] x k+1 = (Φ k ) H b k+1 .

[0028] S10: Calculate the defined transform Ψ in the (k + 1)-th iteration k+1 , and the calculation formula is as follows:

[0029]

[0030] S11: Calculate the tight frame Φ in the (k + 1)-th iteration k+1 , and the calculation formula is as follows:

[0031]

[0032] S12: For the reconstructed multi-component image x k +1, use the Square Root of Sum of Squares (SOS) to calculate the reconstructed single-coil magnitude image, and the calculation formula is as follows:

[0033]

[0034] where represents the single-coil magnitude reconstruction image in the (k + 1)-th iteration.

[0035] S13: Calculate the first intermediate variable in the (k + 1)-th iteration

[0036] S14: Calculate the second intermediate variable in the (k + 1)-th iteration

[0037] S15: Calculate the relative error (RE) between and , and the calculation formula is:

[0038]

[0039] where represents the single-coil magnitude reconstruction image in the k-th iteration.

[0040] S16: Determine whether the algorithm stop criterion is satisfied. If RE < tol is satisfied, go to step S17; otherwise, set k = k + 1 and return to S1;

[0041] S17: Output the final reconstructed single-coil magnitude image

[0042] The beneficial effects of the present invention are as follows: the ESPIRiT model is a parallel magnetic resonance imaging reconstruction technology that utilizes iterative self-consistency of eigenvector operators, and has strong flexibility and compatibility in magnetic resonance image reconstruction. The present invention introduces data-driven adaptive sparse transform learning into the ESPIRiT model, and proposes a parallel MRI reconstruction algorithm based on ESPIRiT multiple groups of sensitivity operators combined with data-driven adaptive sparse transform learning regularization terms, which further enhances the performance of the ESPIRiT model and ensures the quality of MRI reconstructed images. The algorithm transforms the reconstruction problem into a minimization problem containing composite terms, obtains a series of unconstrained problems through the Bregman iteration technique, and transforms the unconstrained problem into a gradient problem and a denoising problem for solution using the operator splitting technique. Theoretical analysis and experimental results show that compared with other algorithms based on the ESPIRiT model (such as TV-ESPIRiT and LpJTV-ESPIRiT), the TL-ESPIRiT algorithm proposed in the present invention has better denoising and restoration effects on images, reduces the staircase artifacts of the reconstructed image, retains the texture details intact, and has a high consistency between the target image and the original image. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 is a flow chart of the method of the present invention;

[0044] Figure 2 To obtain the in vivo human brain slice data of the subject using a 12-channel head coil for full sampling, the 95th slice was taken with a size of 218 × 170 (i.e., dataset 1);

[0045] Figure 3 It is a 2D Poisson disk undersampling mask with a 4x acceleration factor and a 24×24 center fully sampled self-calibration area;

[0046] Figures 4 to 6 The images reconstructed from dataset 1 using the TV-ESPIRiT algorithm, LpJTV-ESPIRiT algorithm, and TL-ESPIRiT algorithm, respectively, undersampled from a two-dimensional Poisson disk sampling mask with a 4x acceleration factor and a 24×24 center fully sampled self-calibration area;

[0047] Figures 7 to 9 They are Figures 4 to 6 The error map corresponding to the reconstructed image, that is, the error map reconstructed by the TV-ESPIRiT algorithm, the LpJTV-ESPIRiT algorithm, and the TL-ESPIRiT algorithm. DETAILED DESCRIPTION

[0048] The technical solution of the present invention is further described in detail below with reference to the accompanying drawings:

[0049] Embodiment 1: The present invention is an efficient reconstruction method proposed based on the ESPIRiT framework.

[0050] 1) ESPIRiT framework:

[0051] Assumptions Represents the multi-component image to be reconstructed, N = N h ×N v Represents the number of pixels of the single coil amplitude image, N v and N h are the number of rows and columns of the single coil amplitude image, x j represents the jth component image to be reconstructed, j=1...J represents the index of the sensitivity map group, for example, x 1 and x J They respectively represent the 1st component image and the Jth component image of the column-vectorized multi-component image x to be reconstructed, J represents the total number of groups of sensitivity maps, and J is 2 in the present invention. represents undersampled multi-coil K-space data, represents the undersampled K-space data of the cth coil, M is the actual number of sampling points of the single coil K-space data, and M<<N, c=1......C represents the index of the receiving coil, for example, y 1 and C They respectively represent the 1st coil data and the Cth coil data of the column-vectorized under-sampled multi-coil K-space data y, and C represents the number of receiving coils used for parallel imaging.

[0052] Then, the ESPIRiT reconstruction model is obtained as:

[0053] y=RFSx (1)

[0054] in, is an undersampling matrix, represents the matrix for selecting sampling point locations from the K-space grid, I C is the C×C identity matrix, represents the Kronecker product. is the Fourier operator, and Respectively represent N h point and N, the two-dimensional Fourier transform matrix of the point. is the sensitivity map group matrix, obtained by the ESPIRiT model, and the j-th group sensitivity map of the c-th coil of S is expressed by S cj Indicates, for example: S 11 The first group of sensitivity diagrams for the first coil, S CJRepresents the Jth group of sensitivity diagrams for the Cth coil.

[0055] In ESPIRiT, the calculation formula for the sensitivity operator is:

[0056]

[0057] Among them, s r yes The eigenvalue of is "1", and r represents the k-space position index variable; Represents the convolution of a semidefinite matrix value at each position r, defined as is a convolution with a matrix-valued kernel, defined as:

[0058]

[0059] Among them, V || It is obtained from the k-space calibration matrix. If the calibration matrix is ​​defined as A, the singular value decomposition of the matrix A is:

[0060] A=U∑V H (4)

[0061] Where Λ is the singular value diagonal matrix of A, U and V are both unitary matrices, whose columns are the left and right singular vectors of the singular values. Decomposing V can obtain the null space span V of A. ⊥ and the row space span V of A || Since the columns of V are the basis of the rows of matrix A, all the self-calibration information is located in V || middle.

[0062] 2) Algorithm derivation:

[0063] Although the ESPIRiT algorithm combined with the LpJTV regularization method can improve the reconstruction quality to a certain extent, there is still much room for improvement. In order to better study the reconstruction performance of ESPIRiT, reduce the staircase artifacts of the reconstructed image, and better preserve the edge contour information. The present invention introduces data-driven adaptive sparse transform learning into ESPIRiT and uses two sets of sensitivity information operators for image reconstruction. The reconstruction problem is then expressed as:

[0064]

[0065] in, Extract the matrix for the image block, which can extract the column-vectorized Image blocks. Represents the sparse transformation matrix of the image block, subject to the constraint W T W=I;||·|| 0Indicates the number of vector or matrix elements with non-zero modulus (i.e., L0 norm). α is the regularization parameter, N P Indicates the number of image blocks into which a single coil image is divided.

[0066] Define the transformation: Easy to get: H Ψ=nI N .make: Then we have: Φ H Φ=I 2N , which shows that Φ is a tight frame, then we can get: Formula (5) can be rewritten as the following optimization problem:

[0067]

[0068] Since Φ is a tight frame, using the tight frame pFISTA technique, the optimization problem (6) is transformed into the following sub-problems:

[0069]

[0070]

[0071] x k+1 =(Φ k ) H b k+1 (9)

[0072] in: L is The Lipschitz constant of the gradient of , λ = α / L represents the parameter of the regularization term. Formula (8) can be further transformed into:

[0073]

[0074] make: Then we have:

[0075]

[0076] make:

[0077]

[0078] Then we get: So formula (10) can be rewritten as:

[0079]

[0080] Z k+1 (B k ) H Perform singular value decomposition to get Z k+1 (B k )H =U∑V H , then the solution of formula (14) is:

[0081] W k+1 =VU H (15)

[0082] In addition, the solution of equation (11) can be obtained by using the hard threshold method:

[0083] b k+1 =H(Φ k z k+1 ,λ) (16)

[0084] Where H(x, θ) is defined as follows:

[0085]

[0086] x represents the input vector and θ represents the threshold.

[0087] In summary, all sub-problems can be solved effectively. A new ESPIRiT parallel MRI reconstruction algorithm with TL regularization term is obtained, which is called TL-ESPIRiT algorithm.

[0088] Reconstruct the multi-component image x k+1 The reconstructed single coil amplitude image is calculated using the square root of the sum of squares (SOS):

[0089]

[0090] Among them, among them, It represents the single coil amplitude reconstruction image of the kth iteration. When the relative error (RE) is lower than the tolerance tol, the algorithm stops. and The RE is defined as:

[0091]

[0092] The specific flow chart is as follows Figure 1 As shown, the steps are as follows:

[0093] S0: Initialization, let k = 0, x 0 =(RFS) H y, t 0 =0, b 0 =H(Φ 0 z 0 ,λ),

[0094] S1: Calculate the intermediate variable z for the k+1th iteration k+1 , the calculation formula is as follows (7);

[0095] S2: Initialize j=1;

[0096] S3: Update B k = unvec(b k ), the calculation formula is as follows (12);

[0097] S4: Calculate the intermediate variable Z for the k+1th time k+1 , the calculation formula is as follows (13);

[0098] S5: Z k+1 (B k ) H Perform singular value decomposition and get Z k+1 (B k ) H =U∑V H ;

[0099] S6: Transform update, calculate the sparse transformation W of the k+1th iteration k+1 , the calculation formula is as follows (15);

[0100] S7: Hard threshold denoising, calculate the auxiliary variable b for the k+1th iteration k+1 , the calculation formula is as follows (16);

[0101] S8: Determine whether all image components are completed. If j=J, go to step S9; otherwise, set j=j+1 and return to S3.

[0102] S9: Calculate the image x to be reconstructed at the k+1th iteration k+1 , the calculation formula is as follows (9);

[0103] S10: Calculate the definition transformation Ψ for the k+1th iteration k+1 , the calculation formula is as follows

[0104] S11: Calculate the tight frame Φ for the k+1th iteration k+1 , the calculation formula is as follows

[0105] S12: Calculate the reconstructed single coil amplitude image, the calculation formula is as shown in (18);

[0106] S13: Calculate the first intermediate variable t of the k+1th iteration k+1 , the calculation formula is as follows

[0107] S14: Calculate the second intermediate variable of the k+1th iteration The calculation formula is as follows

[0108] S15: Calculation and The relative error (RE) between them is calculated as shown in (19);

[0109] S16: Determine whether the algorithm stopping criteria are met. If RE<tol is met, go to step S17; otherwise, set k=k+1 and return to S1;

[0110] S17: Output the final reconstructed single coil amplitude image Experimental results:

[0111] In the following experiments, in order to verify the reconstruction performance of the TL-ESPIRiT algorithm, the proposed TL-ESPIRiT algorithm is compared with the TV-ESPIRiT and LpJTV-ESPIRiT methods. All algorithms can obtain reconstructed images, and all experiments are implemented on Matlab (MathWorks, Natick, MA). All experiments are performed on a laptop computer with an Intel Core i5-4210U @ 2.4 GHz processor, 8 GB memory, and Windows 10 operating system (64-bit).

[0112] In order to compare the performance of each algorithm, the present invention uses a human brain slice data for simulation experiments, named dataset 1. Figure 2 As shown in Figure 1, dataset 1 is a 12-channel head coil fully sampled to obtain the in vivo human brain slices of the subject, taking the 95th slice with a size of 218×170. A two-dimensional Poisson disk undersampling mask with an acceleration factor of AF=4 and a 24×24 center fully sampled self-calibration area is used, as shown in Figure 3 As shown, artificial sampling is performed on the full sampling data to generate test data.

[0113] The parameters of all algorithms are adjusted for the optimal signal-to-noise ratio (SNR). All experiments use SNR to quantitatively judge the reconstruction quality. The SNR value is calculated in the region of interest. The higher the SNR, the better the reconstruction effect.

[0114] SNR is defined as follows:

[0115]

[0116] Among them, MSE represents the reconstructed image and the reference image x, Var represents the variance of x.

[0117] First, the three reconstruction algorithms under dataset 1 are visually compared. The present invention selects the reconstructed images when the acceleration factor is 4 for comparison. Figures 4 to 6 shown. Figures 4 to 6 The reconstructed images of the algorithms TV-ESPIRiT, LpJTV-ESPIRiT and TL-ESPIRiT are shown in turn. Figure 4 and Figure 5 There are blurring artifacts and some details are not reconstructed. Figure 5 Although there are a small amount of artifacts, the artifact removal effect is significantly improved compared with the previous two, and the reconstructed image is more consistent with the original image.

[0118] In order to further illustrate that the new algorithm TL-ESPIRiT proposed in this invention is superior to the algorithms TV-ESPIRiT and LpJTV-ESPIRiT in terms of reconstruction performance, Figures 7 to 9 The error graphs of dataset 1 under the three algorithms are shown (the whiter the error graph, the greater the error). It can be seen from the error graph that the TL-ESPIRiT algorithm has a smaller reconstruction error, and the reconstruction performance of the algorithm proposed in the present invention is better. In addition, the experiment quantitatively evaluates the performance of the three reconstruction algorithms through the SNR value. The algorithm TL-ESPIRiT proposed in the present invention significantly improves the SNR value of the reconstructed image.

[0119] In summary, experiments were conducted on different reconstruction algorithms for the selected dataset dataset 1. The TL-ESPIRiT algorithm performed well in both evaluation indicators and visual aspects, and the reconstructed image quality was better than that of the TV-ESPIRiT and LpJTV-ESPIRiT algorithms.

[0120] The specific implementation modes of the present invention are analyzed and described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above implementation modes, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. An improved ESPIRiT reconstruction method based on sparse transform learning, Features: The following steps are involved: S0: Initialization, let k = 0, x 0 =(RFS) H y, t 0 =0, b 0 =H(Φ 0 z 0 ,λ), Where k is the loop variable, indicating the kth iteration of the data, and the superscript " 0 " indicates the initial value, Represents the multi-component image to be reconstructed, N = N h ×N v Represents the number of pixels of the single coil amplitude image, N v and N h are the number of rows and columns of the single coil amplitude image, x j represents the jth component image to be reconstructed, j=1...J represents the index of the sensitivity map group, x 1 and x J They represent the first component image and the Jth component image of the column vectorized multi-component image x to be reconstructed, J represents the total number of sensitivity maps, x 0 represents the initial value of x, represents undersampled multi-coil K-space data, represents the undersampled K-space data of the cth coil, M is the actual number of sampling points of the single coil K-space data, and M<<N, c=1......C represents the index of the receiving coil, y 1 and C They represent the first coil data and the Cth coil data of the column-vectorized undersampled multi-coil K-space data y, respectively. C represents the number of receiving coils used for parallel imaging. The superscript " T " indicates matrix transpose, superscript " H " represents the complex conjugate transpose of the matrix, is an undersampling matrix, where represents the matrix for selecting sampling point locations from the K-space grid, I C is the C×C identity matrix, represents the Kronecker product, is the Fourier operator, and Respectively represent N h Point and N v The two-dimensional Fourier transform matrix of a point, is the sensitivity map group matrix, obtained by the ESPIRiT model, and the j-th group sensitivity map of the c-th coil of S is expressed by S cj Indicates that S 11 The first group of sensitivity diagrams for the first coil, S CJ represents the J-th group sensitivity diagram of the C-th coil, and represents the intermediate variable, t 0 and are t and The corresponding initial value; P i (·): It represents the linear operator extracted from the i-th image block, which can extract J images of size The extracted J image blocks are converted into The matrix, P 1 (·) indicates the first matrix extraction linear operator, Indicates the Nth P The matrix extracts the linear operator, express The sparse transformation matrix of the image patch, W 0 represents the initial value of W, express The discrete cosine transform matrix of the point, Ψ is the defined transformation, and Ψ 0 is the initial value of Ψ; Φ is a tight frame, and Φ H Φ=I 2nN , Φ 0 is the initial value of Φ; H(x, θ) is the hard threshold function, x represents the input matrix, θ represents the threshold, λ is the parameter, and λ>0, z is the intermediate variable, and all z are concatenated into a matrix in order And Z represents 2N image blocks extracted from z and spliced ​​into an n×2N matrix, P 1 z 1 and P N z 2 They represent the first image block and the 2Nth image block extracted from the multi-component image z by the intermediate variables, respectively. It means Z=(P 1 z 1 … P N z 1 P 1 z 2 … P N z 2 ), the initial value of auxiliary variables is the column-vectorized encoding matrix, b 0 represents the initial value of b; S1: Calculate the intermediate variable z for the k+1th iteration k+1 , the calculation formula is as follows: in, Represents the intermediate variable of the kth iteration L is The Lipschitz constant of the gradient; S2: Initialize j=1; S3: Update B k = unvec(b k ), the calculation formula is as follows: S4: Calculate the k+1th intermediate variable, Z k+1 : S5: Z k+1 (B k ) H Perform singular value decomposition and get Z k+1 (B k ) H =U∑V H ; Among them, B k is the kth iteration of the auxiliary variable B, Z k+1 is the k+1th iteration of the intermediate variable Z, U, ∑, V are the singular value decomposition factors; S6: Transform update, calculate the sparse transformation W of the k+1th iteration k+1 , the calculation formula is as follows: W k+1 =VU H ; S7: Hard threshold denoising, calculate the auxiliary variable b for the k+1th iteration k+1 , the calculation formula is as follows: b k+1 =H(Φ k z k ,l); S8: Determine whether all image components are completed. If j=J, go to step S9; otherwise, set j=j+1 and return to S3; S9: Calculate the image x to be reconstructed at the k+1th iteration k+1 , the calculation formula is as follows: x k+1 =(Φ k ) H b k+1 ; S10: Calculate the definition transformation Ψ for the k+1th iteration k+1 , the calculation formula is as follows: S11: Calculate the tight frame Φ for the k+1th iteration k+1 , the calculation formula is as follows: S12: Reconstruct the multi-component image x k+1 The square root of the sum of squares SOS is used to calculate the reconstructed single-coil amplitude image, and the calculation formula is as follows: in, represents the single coil amplitude reconstructed image of the k+1th iteration; S13: Calculate the first intermediate variable of the k+1th iteration S14: Calculate the second intermediate variable of the k+1th iteration S15: Calculation and The relative error RE between them is calculated as follows: in, represents the single coil amplitude reconstructed image of the kth iteration; S16: Determine whether the algorithm stopping criteria are met. If RE<tol is met, go to step S17; otherwise, set k=k+1 and return to S1; S17: Output the final reconstructed single coil amplitude image

Images