CN106019189A - Self consistency based parallel magnetic resonance imaging quick reconstructing method - Google Patents

Abstract

The invention discloses a self consistency based parallel magnetic resonance imaging quick reconstructing method. Based on an SPIRiT frame, aiming at a JTV and JL1 composite regular term contained parallel imaging reconstruction problem, the invention proposes a quick reconstructing method employing Descartes sampling. In the prior art, scanning time can be reduced distinctively by adopting the magnetic resonance imaging technology; however, the reconstructing time of a prior reconstructing algorithm is long. The method provided by the invention utilizing a split Bregman technology for dissolving a prior problem into a plurality of sub problems that are easy to solve. Experiment simulation results show that the new algorithm is greatly improved in convergence speed under the premise of ensuring reconstruction quality compared with a prior NLCG algorithm.

Description

A kind of based on from conforming parallel MR imaging method for fast reconstruction

Technical field

The present invention relates to a kind of based on from conforming parallel MR imaging method for fast reconstruction.

Background technology

Nuclear magnetic resonance (Magnetic Resonance Imaging, MRI) can provide good human body soft tissue pair Ratio degree, and there is no radioactivity, therefore it is increasingly becoming indispensable imaging tool in current clinical medicine.But, by physics With the restriction of physiologic factor, the speed of collecting magnetic resonance signal is the slowest.Parallel imaging is a kind of conventional skill improving picking rate Art.The proposition of SMASH and SENSE, indicates that parallel imaging becomes practicable technology.Parallel imaging uses to be believed magnetic resonance Number multiple coils that sensitivity is different collecting magnetic resonance signal simultaneously.The use of sensitivity information decreases the data for reconstruct Number, thus improve the speed of collection.

More than ten years in past, create a variety of parallel imaging reconstruction algorithm.The difference of these algorithms is: 1) explicit Or implicit expression uses sensitivity information；2) recover single image or by-line circle recovers multiple image.SMASH, SENSE, SPACE-RIP, kSPA belong to the method using explicit sensitivity information to recover single image.PILS, GRAPPA, SPIRiT belong to Use the method that implicit expression sensitivity information carries out by-line circle recovery.And AUTO-SMASH uses implicit expression sensitivity information to recover single Image, PARS then uses explicit sensitivity information to carry out by-line circle recovery.

SENSE is to apply most common in various algorithm, and it is easily integrated the priori of image.In recent years, there is research Compressed sensing is combined by personnel with SENSE, increases regular terms, achieves preferable quality reconstruction.When sensitivity is known, SENSE Optimal solution can be obtained.But accurately estimate that the sensitivity of coil is the most extremely difficult.Therefore, it is not required to explicitly use coil sensitivity The method for self-calibrating of information has just embodied advantage.

GRAPPA is the self calibration reconstructing method of a kind of widely used by-line circle, sensitive owing to need not explicit use Degree information, thus avoid thorny sensitivity and estimate.Lustig etc. propose one more effectively on the basis of GRAPPA The theoretical frame SPIRiT of parallel imaging technique.Similar with SENSE, reconstruction is changed into one against asking by SPIRiT Topic solves, it is easy to the priori of integrated image, and can reconstruct the data of any k-space sampling.Murphy and Vasanawala etc. propose L1-SPIRiT model, introduce Joint L1 regular terms, and use POCS and NLCG scheduling algorithm to ask Solve.But the reconstruction quality of POCS is not good enough, and the reconstitution time of NLCG algorithm is longer.

Summary of the invention

It is an object of the invention to overcome the deficiencies in the prior art, it is provided that one becomes based on from conforming parallel MR As method for fast reconstruction.Parallel MR imaging technology can significantly reduce sweep time, but the reconstruct of prior art is calculated Method reconstitution time is longer.The present invention is based on SPIRiT framework, for the weight of the parallel imaging being combined regular terms containing JTV and JL1 Structure problem, it is proposed that the fast Reconstruction Algorithms of a kind of Descartes sampling.This algorithm utilizes division Bregman technology (split Bregman) it is multiple subproblems being prone to and solving by former PROBLEM DECOMPOSITION.Simulation results shows, with existing NLCG algorithm Comparing, new algorithm is in the case of ensureing reconstruction quality, and convergence rate has greatly improved.

It is an object of the invention to be achieved through the following technical solutions: a kind of become based on from conforming parallel MR As method for fast reconstruction, it comprises the following steps:

S0: initialize, make x⁰=0,w⁰=0, d⁰=0, z⁰=0, k=0；

In formula,For multi-coil image variables, each it is classified as a coil image by row storehouseN =m × n, C are multi-coil number, m and n is respectively line number and the columns of single coil two dimensional image；K is cyclic variable；b_w、b_d And b_zFor dual variable, w, d and z are aleatory variable, and w=x, z=Ψ x, Ψ are the wavelet transformation of by-line circle, Wherein D_nAnd D_mRepresenting the circular matrix of n × n and m × m respectively, described circular matrix structure is as follows:

S1: calculate w respectively^k+1、d^k+1And z^k+1Value, wherein:

$w^{k+1}={\mathrm{\Γ}}^{-1}\left(\mathrm{\μ}\right(x^k+b_w^k\left)\right);$

In formula, matrix Γ=λ (G-I)^T(G-I)+μ I has diagonal blocks structure, simply reorders this matrix all Block diagonal arrangement can be produced, be made up of C × C block；μ is punishment parameter, and λ is punishment parameter；

$z^{k+1}=shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1);$

In formula, μ₁Representing punishment parameter, shrinkJ () represents the one-dimensional contraction operator of associating, and computing formula is as follows:

$shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1)=max\left(\right|{\mathrm{\Ψx}}^k+b_z^k|-1/{\mathrm{\μ}}_1,0)\·\frac{{\mathrm{\Ψx}}^k+b_z^k}{s_1};$

In formula,Ψ represents the wavelet transformation of by-line circle；Repmat is in matlab Function, function be by first variable according to the dimension of second parameter carry out replicate expand, c represents loop index；

$d^{k+1}=shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2);$

In formula, μ₂Representing punishment parameter, shrink2J () represents associating two-dimensional contraction operator, and computing formula is as follows:

$shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2)=max\left(\right|{\mathrm{Dx}}^k+b_d^k|-1/{\mathrm{\μ}}_2,0)\·\frac{{\mathrm{Dx}}^k+b_d^k}{s_2};$

In formula,

S2: in the case of the value of cyclic variable c is respectively 1～C, calculates

In formula,For Fourier transform matrix,Represent coefficient of frequency select son, its every a line all from Extracting out in the unit matrix of N × N, M is that the frequency of lack sampling is counted, and M/N is lack sampling rate； For adopting The coefficient of frequency that collection arrives, its every string is the lack sampling coefficient of frequency of a coil image, α₁And α₂For regularization parameter；

S3: the most right WithValue be updated, wherein:

$b_w^{k+1}=b_w^k+(x^{k+1}-w^{k+1});$

$b_z^{k+1}=b_z^k+({\mathrm{\Ψx}}^{k+1}-z^{k+1});$

$b_d^{k+1}=b_d^k+({\mathrm{Dx}}^{k+1}-d^{k+1});$

S4: judge that k value, whether more than maximum cycle K, returns step if not after k value then adds 1 operation Rapid S1, otherwise enters step S5；

S5: the x that output finally gives:

X=x^k+1；

S6: the view data for each coil at each loop restructuring uses SRSOS method to each coil image Combine, obtain final single width reconstruct image:

$SRSOS\left(x\right)=\sqrt{{\Σ}_c{\left|x_c\right|}^2}.$

The invention has the beneficial effects as follows: the invention discloses that one quickly weighs based on from conforming parallel MR imaging Structure method.Parallel MR imaging technology can significantly reduce sweep time, but existing restructing algorithm reconstitution time is relatively Long.The present invention is based on SPIRiT framework, for the reconstruction of the parallel imaging being combined regular terms containing JTV and JL1, it is proposed that A kind of fast Reconstruction Algorithms of Descartes sampling.This algorithm utilizes division Bregman technology (split Bregman) by former problem It is decomposed into multiple subproblem being prone to and solving.Simulation results shows, compared with existing NLCG algorithm, new algorithm is ensureing In the case of reconstruction quality, convergence rate has greatly improved.

Accompanying drawing explanation

Fig. 1 is the inventive method flow chart；

Fig. 2 is the fully sampled high-resolution brain imaging data (i.e. data1) of the use SPGR retrieval of T1 weighting Figure；

Fig. 3 is fully sampled axial Brian Imaging data (i.e. data2) figure using 2 dimensions to rotate back to city retrieval；

Fig. 4 is the Poisson sub-sampling schematic diagram that speed-up ratio is about 6；

Fig. 5 is the sub-sampling rate present invention and NLCG Performance comparision when being the employing data1 cycle tests of 6；

Fig. 6 is the sub-sampling rate present invention and NLCG Performance comparision when being the employing data2 cycle tests of 6；

Fig. 7 is cycle tests error image of NLCG algorithm reconstruct image and original image when being data1；

Fig. 8 is cycle tests error image of the inventive method reconstruct image and original image when being data1；

Fig. 9 is cycle tests error image of NLCG algorithm reconstruct image and original image when being data2；

Figure 10 is cycle tests error image of the inventive method reconstruct image and original image when being data2.

Detailed description of the invention

Technical scheme is described in further detail below in conjunction with the accompanying drawings:

The present invention is a kind of efficient reconstructing method proposed based on SPIRiT framework.

SPIRiT be one based on self-alignment parallel imaging reconstruction method, lack when the method by-line circle interpolation goes out sub-sampling The Frequency point lost.In SPIRiT, an interpolation core g_ijBy data (the commonly referred to self-correcting definite message or answer fully sampled to frequency domain center Number) carry out calibration and obtain.If with x_iRepresent the whole frequency domain data of i-th coil, then conformance criteria based on calibration is writeable Become:

$x_i=\underset{j=1}{\overset{C}{\Σ}}{g}_{ij}*x_j---\left(1\right)$

In formula, C represents coil number, and " * " represents convolution operation.Convolution kernel g_ijIt is referred to as SPIRiT core.

All the conformance criteria of coils can be simplified to matrix form:

X=Gx (2)

In formula, x is the view data of whole coil.G is interpolating operations, it is possible to as the convolution of frequency domain or image area Product realizes.

Certainly, except concordance based on calibration, reconstruct also needs to the data with former sampling and keeps concordance, is expressed as square Formation formula:

WhereinFor multi-coil image variables, each it is classified as a coil image by row storehouseN= M × n, m and n are respectively line number and the columns of single coil two dimensional image.For Fourier transform matrix, Representing that coefficient of frequency selects son, its every a line is all extracted out from the unit matrix of N × N, and M is that the frequency of lack sampling is counted, and M/N is Lack sampling rate； For the coefficient of frequency collected, its every string is the lack sampling of a coil image Coefficient of frequency.

Owing to multi-coil image divides the similarity of transform domain at wavelet transformed domain and first difference, therefore we draw simultaneously Enter Joint L1 (JL1) regular terms and Joint Total Variation (JTV) regular terms.According to formula (2) and the two of formula (3) Individual restriction, the reconstruction of parallel imaging can change into following optimization problem:

Wherein α₁And α₂For regularization parameter, Ψ represents the wavelet transformation of by-line circle, ||x||_JTVIt is defined as the following formula:

${||x||}_{JTV}=\underset{r=1}{\overset{N}{\Σ}}\sqrt{\underset{c=1}{\overset{C}{\Σ}}\[{\left|{\left(D_{h}x\right)}_{rc}\right|}^2+{\left|{\left(D_{v}x\right)}_{rc}\right|}^2\]}---\left(5\right)$

Wherein, D_hX and D_vX respectively carries out the first difference of line direction and column direction and divides conversion x, D_mRepresenting the circular matrix of m × m, structure is as follows:

In order to derive conveniently, the first difference in row, column direction divide conversion merge and writes: I.e.Therefore, the first difference of x divides conversion can be shown as D by simple table_hvx.Clearly have:

Derive for convenience, by | | x | |_JTVWrite as following form:

${||x||}_{JTV}={||D_{hv}x||}_{g,1}=\underset{r=1}{\overset{N}{\Σ}}\sqrt{\underset{c=1}{\overset{C}{\Σ}}\[{\left|{\left(D_{h}x\right)}_{rc}\right|}^2+{\left|{\left(D_{v}x\right)}_{rc}\right|}^2\]}---\left(7\right)$

Therefore, the optimization problem shown in formula (4) can be write as:

In order to solve the optimization problem shown in formula (8), apply quadratic penalty function technology, then can obtain:

Wherein λ is punishment parameter.

In order to solve the optimization problem shown in formula (9), introducing aleatory variable w=x, z=Ψ x, d=Dx, d have following shape Formula:

$d=\left(\begin{array}{c}{d}_h\\ d_v\end{array}\right)=\left(\begin{array}{c}{D}_{h}x\\ D_{v}x\end{array}\right)=D_{hv}x=\left(\begin{array}{ccc}{d}_{h11}& \mathrm{...}& d_{h1C}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ d_{hN1}& \mathrm{...}& d_{hNC}\\ d_{v11}& \mathrm{...}& d_{v1C}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}\\ d_{vN1}& \mathrm{...}& d_{vNC}\end{array}\right)---\left(10\right)$

$d={(D_h^T,D_v^T)}^T={({\left(D_{h}x\right)}^T,{\left(D_{v}x\right)}^T)}^T=Dx.$

Utilizing division Bregman technology, former problem formula (9) can pass through following iterative problem solving:

$b_w^{k+1}=b_w^k+(x^{k+1}-w^{k+1})---\left(12\right)$

$b_z^{k+1}=b_z^k+({\mathrm{\Ψx}}^{k+1}-z^{k+1})---\left(13\right)$

$b_d^{k+1}=b_d^k+({\mathrm{Dx}}^{k+1}-d^{k+1})---\left(14\right)$

In formula (11), μ, μ₁And μ₂It is punishment parameter.Formula (11) is solved x, w, z, d respectively, then obtains following subproblem Solve:

$w^{k+1}=\arg \underset{x}{\min }\frac{\mathrm{\μ}}{2}{||w-x^k-b_w^k||}_2^2+\frac{\mathrm{\λ}}{2}{||(G-I)w||}_2^2---\left(15\right)$

$z^{k+1}=\arg \underset{z}{\min }\frac{{\mathrm{\μ}}_1}{2}{||z-{\mathrm{\Ψx}}^k-b_z^k||}_2^2+{||z||}_{2,1}---\left(16\right)$

$d^{k+1}=\arg \underset{d}{min}\frac{{\mathrm{\μ}}_2}{2}{||d-{\mathrm{Dx}}^k-b_d^k||}_2^2+{||d||}_{g,1}---\left(17\right)$

The optimal condition of formula (15) isAnd matrix

Γ=λ (G-I)^T(G-I)+μ I has diagonal blocks (diagonal-block) structure, carries out this matrix simply Reorder and all can produce block diagonal angle (block-diagonal) structure, be made up of C × C block.The inverse renewal occurring in w of Γ is expressed In formula:

$w^{k+1}={\mathrm{\Γ}}^{-1}\left(\mathrm{\μ}\right(x^k+b_w^k\left)\right)---\left(19\right)$

Owing to C is smaller, therefore directly each C × C block in Γ can be inverted and try to achieve Γ^-1。

In formula (16), the solution of z is given by:

$z^{k+1}=shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1)---\left(20\right)$

In formula, μ₁Representing punishment parameter, shrinkJ () represents the one-dimensional contraction operator of associating (combining one-dimensional Soft thresholding), Computing formula is as follows:

$shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1)=max\left(\right|{\mathrm{\Ψx}}^k+b_z^k|-1/{\mathrm{\μ}}_1,0)\·\frac{{\mathrm{\Ψx}}^k+b_z^k}{s_1}---\left(21\right)$

In formula,Repmat is the function in matlab, and function is by first change Measure to carry out replicating according to the dimension of second parameter and expand.

In formula (17), the solution of d is given by:

$d^{k+1}=shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2)---\left(22\right)$

In formula, μ₂Representing punishment parameter, shrink2J () represents associating two-dimensional contraction operator (associating two dimension Soft thresholding), Computing formula is as follows:

$shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2)=max\left(\right|{\mathrm{Dx}}^k+b_d^k|-1/{\mathrm{\μ}}_2,0)\·\frac{{\mathrm{Dx}}^k+b_d^k}{s_2}---\left(23\right)$

In formula,

Owing in formula (18), each coil variable is independent, therefore single coil can be separated into and calculate, it may be assumed that

According to optimal condition, the solution of formula (24) can be obtained by following formula:

Owing to Ψ is orthogonal wavelet, therefore there is Ψ^TΨ=I；Can be by FFT diagonalization, therefore this formula Can pass through FFT rapid solving:

Since all subproblems can Efficient Solution, for solving of problem (9), the present invention obtains a kind of efficient weight Structure algorithm quickly divides Bregman algorithm (Split based on from conforming containing the parallel MR imaging combining full variation Bregman for SPIRiT based Parallel MR Imaging with JTV and JL1regularization, SB4SpMRI), the idiographic flow of SB4SpMRI algorithm is as shown in Figure 1.In this flow process, C is coil sum, and c is for follow accordingly Ring variable (the most previously described coil index)；K is total iterations, and k is corresponding cyclic variable (i.e. kth time iteration).

In algorithm 1, the 2nd step is loop control condition, when iteration reaches maximum iteration time or when reconstruct image is with front When the relative error of the image that secondary iterative reconstruction goes out is less than certain value, stop circulation.In algorithm 1, what each loop restructuring went out is The view data of each coil, therefore also needs to use SRSOS (Square Root Sum Of Squares) method to each line Loop graph picture is combined, just available final single width reconstruct image.This last handling process can represent with SRSOS (x), calculates Shown in method such as formula (27):

$SRSOS\left(x\right)=\sqrt{{\Σ}_c{\left|x_c\right|}^2}---\left(27\right)$

Algorithm described herein use matlab programmed environment (Version R2008b, the MathWorks Inc., Natick, MA) realize, and apply in the parallel MR image reconstruction of sub-sampling.All experiments are being configured to Intel Core Perform on CPU, the 4GB internal memory of i53230M@2.6GHz and the notebook of Windows 764bit operating system.

Material used in experiment has two groups, is all real multi-coil concurrent MR imaging data, as shown in Figures 2 and 3:

First group is the imaging data of fully sampled high-resolution brain of use SPGR retrieval of T1 weighting (data1).Using has the GE Signa-Excite 1.5T scanning system of 8 channel reception coils to scan.Sweep parameter is: TE and TR is respectively 8ms and 17.6ms, flip and is equal to 20 °, BW=6.94kHz.FOV is 20 × 20 × 20cm, and imaging array is big Little is 200 × 200 × 200.

Second group is fully sampled axial Brian Imaging data (data2).This imaging data uses 2 dimension spin-echo sequences (2D spin echo sequence) is having the GE 3T scanning device (GE Healthcare, Waukesha, WI) of 8 coils Up-sampling obtains, TE and TR is respectively 11ms and 700ms, and matrix size is 256 × 256, and FOV is 220 × 220mm.

The present invention uses the Poisson sub-sampling pattern having been widely used that fully sampled data are carried out sub-sampling, thus obtains not With the sub-sampling data accelerating multiple.In sub-sampling pattern, for the fully sampled signal in the self-alignment center of SPIRiT core a size of 24×24.Fig. 4 is that speed-up ratio is about the schematic diagram of the sample pattern of 6 (the fully sampled signal in the center that has contemplated that, matrix size is 256 ×256)。

Being used herein as SNR and carry out the quality of quantitative assessment reconstruct image, SNR is defined as:

$SNR=10{\log }_{10}\left(\frac{Var}{MSE}\right);---\left(34\right)$

In formula, MSE represents the mean square error between reference picture and reconstruct image, and Var represents the variance of reference picture.

In order to test the performance of the new algorithm of the present invention, we are by new algorithm (SB4SpMRI) and classic algorithm non-thread The algorithm (NLCG) of property conjugate gradient compares, and regular terms all uses JTV and JL1 to be combined regular terms.SPIRiT core in experiment Size is all provided with being set to 5 × 5.The regularization parameter of each algorithm is carried out tuning, so that algorithm reaches optimal performance.

Fig. 5 and Fig. 6 gives the Performance comparision of new algorithm SB4SpMRI Yu NLCG when sub-sampling rate is 6.When certain algorithm When the SNR change of reconstruct image tends towards stability, it is considered as this algorithm and has restrained.Fig. 5 and Fig. 6 respectively uses data sequence Experimental result during data1 and data2.From fig. 5, it can be seen that new algorithm SB4SpMRI can reach suitable with NLCG algorithm SNR, and convergence rate is greatly faster than NLCG.And for Fig. 6, it is also possible to obtain the conclusion being similar to.

Fig. 7 and Fig. 8 sets forth cycle tests when being data1, NLCG and SB4SpMRI algorithm reconstruct image is with original The error image of image.Fig. 9 and Figure 10 sets forth cycle tests when being data2, NLCG and SB4SpMRI algorithm reconstruct figure As the error image with original image.The reconstruct image that it is pointed out that the algorithm compared here is all through abundant iteration Obtain afterwards.From Fig. 7, Fig. 8, Fig. 9 and Figure 10 it can be seen that the reconstructed error of SB4SpMRI algorithm is suitable with NLCG.

Claims (1)

1. one kind based on from conforming parallel MR imaging method for fast reconstruction, it is characterised in that: it comprises the following steps:

S0: initialize, make x⁰=0,w⁰=0, d⁰=0, z⁰=0, k=0；

In formula,For multi-coil image variables, each it is classified as a coil image by row storehouseN=m × N, C are multi-coil number, m and n is respectively line number and the columns of single coil two dimensional image；K is cyclic variable；b_w、b_dAnd b_zFor Dual variable, w, d and z are aleatory variable, and w=x, z=Ψ x, Ψ are the wavelet transformation of by-line circle,Wherein D_nAnd D_mRepresenting the circular matrix of n × n and m × m respectively, described circular matrix structure is as follows:

S1: calculate w respectively^k+1、d^k+1And z^k+1Value, wherein:

$w^{k+1}={\mathrm{\Γ}}^{-1}\left(\mathrm{\μ}\right(x^k+b_w^k\left)\right);$

In formula, matrix Γ=λ (G-I)^T(G-I)+μ I has diagonal blocks structure, simply reorders this matrix and all can produce Raw block diagonal arrangement, is made up of C × C block；μ is punishment parameter, and λ is punishment parameter；

$z^{k+1}=shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1);$

In formula, μ₁Representing punishment parameter, shrinkJ () represents the one-dimensional contraction operator of associating, and computing formula is as follows:

$shrinkJ({\mathrm{\Ψx}}^k+b_z^k,1/{\mathrm{\μ}}_1)=max\left(\right|{\mathrm{\Ψx}}^k+b_z^k|-1/{\mathrm{\μ}}_1,0)\·\frac{{\mathrm{\Ψx}}^k+b_z^k}{s_1};$

In formula,Ψ represents the wavelet transformation of by-line circle；Repmat is the letter in matlab Number, function is to carry out replicating according to the dimension of second parameter by first variable expanding, and c represents loop index；

$d^{k+1}=shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2);$

In formula, μ₂Representing punishment parameter, shrink2J () represents associating two-dimensional contraction operator, and computing formula is as follows:

$shrink2J({\mathrm{Dx}}^k+b_d^k,1/{\mathrm{\μ}}_2)=max\left(\right|{\mathrm{Dx}}^k+b_d^k|-1/{\mathrm{\μ}}_2,0)\·\frac{{\mathrm{Dx}}^k+b_d^k}{s_2};$

In formula,

S2: in the case of the value of cyclic variable c is respectively 1～C, calculates x_c ^k+1:

In formula,For Fourier transform matrix,Representing that coefficient of frequency selects son, its every a line is all from N × N Unit matrix in extract out, M is that the frequency of lack sampling is counted, and M/N is lack sampling rate；For collecting Coefficient of frequency, its every string is the lack sampling coefficient of frequency of a coil image, α₁And α₂For regularization parameter；

S3: the most rightWithValue be updated, wherein:

$b_w^{k+1}=b_w^k+(x^{k+1}-w^{k+1});$

$b_z^{k+1}=b_z^k+({\mathrm{\Ψx}}^{k+1}-z^{k+1});$

$b_d^{k+1}=b_d^k+({\mathrm{Dx}}^{k+1}-d^{k+1});$

S4: judge that k value, whether more than maximum cycle K, returns step S1 if not after k value then adds 1 operation, Otherwise enter step S5；

S5: the x that output finally gives:

X=x^k+1；

S6: the view data for each coil at each loop restructuring uses SRSOS method to carry out each coil image Associating, obtains final single width reconstruct image:

$SRSOS\left(x\right)=\sqrt{{\mathrm{\Σ}}_c|x_c{|}^2}.$