CN105022010A - Parallel magnetic resonance image reconstruction method based on regularization iteration - Google Patents

Abstract

The invention discloses a parallel magnetic resonance image reconstruction method based on regularization iteration. The method comprises the following steps: firstly, all K spatial data are subjected to simulation sampling, multi-channel parallel coil K spatial center self-calibration ACS row data and undersampled K spatial data are obtained, and spatial sensitivity distribution of each coil and an aliasing image of each coil are calculated respectively; secondly, a reconstruction image is subjected to equation description combined with a least square method theory, based on image domain method during the calculation process of the aliasing images, a regularization method is introduced, an optimal regularization parameter is calculated and obtained by utilization of an L-curve method, and a matrix equation after regularization is obtained; thirdly, the matrix equation after regularization is subjected to iteration reconstruction by utilization of a conjugate gradient iterative method, and a reconstruction image of each coil is obtained. The method can reduce interference of noise to a reconstruction result further, and has characteristics of high signal to noise ratio, few errors, good imaging effects and the like.

Description

Based on the parallel MR image rebuilding method of regularization iteration

Technical field

The present invention relates to mr imaging technique, be specifically related to a kind of parallel MR image rebuilding method based on regularization iteration.

Background technology

Magnetic resonance imaging (magnetic resonance image, MRI) because of tool radiationless, resolution is high, multi-faceted, and the advantages such as multiparameter are used widely at present clinically.The main weak point of conventional magnetic resonance imaging can be summed up as 2 points: one is that required sweep time is longer, and movement of patient produces artifact, affects clinical diagnosis, and is not suitable for patient of special circumstances; Two is to the such as locomotive organ such as heart, abdominal cavity difficult in imaging.

Parallel MR imaging (pMRI) technology is a great technological breakthrough, and it with step collecting magnetic resonance signal, reduces phase encoding step number via multiple receiving coil, in the reduction sampling time, image taking speed is improved greatly.Adopt pMRI can improve image taking speed and image resolution ratio, but this sacrifices that to rebuild the signal to noise ratio (S/N ratio) of image be cost.The subject matter of parallel imaging technique is the reduction of its signal to noise ratio (S/N ratio), and this comes from the minimizing of hits on the one hand; Come from the pathosis of system matrix on the other hand, noise in inversion process is exaggerated, and rebuild signal noise ratio (snr) of image and decline, quality is not ideal enough.

Summary of the invention

For solving the problems of the technologies described above, the present invention proposes a kind of parallel MR image rebuilding method based on regularization iteration, the method is based on Least Square Theory, regularization method is introduced in the process calculated, utilize L-curve method to solve optimum regularization parameter, carry out image reconstruction in conjunction with conjugate gradient iterative procedure.While improving signal to noise ratio (S/N ratio), strengthen the robustness of the method, effectively solve the problems of the technologies described above, reduce noise to large extent to the impact of reconstructed results, improve signal noise ratio (snr) of image.

For achieving the above object, technical scheme of the present invention is as follows:

Based on the parallel MR image rebuilding method of regularization iteration, method comprises the following steps:

1. pair full K space data carries out analog sampling, obtains the K space data of the capable data of multi-channel parallel coil K space center self-calibration ACS and lack sampling, calculates each coil space sensitivity profile and each coil aliased image respectively;

2. in the computation process that coil aliased image is launched, in conjunction with least squares theory, based on image area method, equation description is carried out to reconstruction image, introduce regularization method, utilize L-curve method to calculate optimum regularization parameter, obtain the matrix equation after regularization;

3. utilize conjugate gradient iterative procedure to carry out iterative approximation to the matrix equation after regularization, obtain the reconstruction image of each coil.

A kind of parallel MR image rebuilding method based on regularization iteration of the present invention, based on Least Square Theory, regularization method is introduced in the process calculated, ask with L-curve method and calculate best regularization parameter, image reconstruction is carried out in conjunction with conjugate gradient iterative procedure, inhibit the impact that the abnormal data of coil causes well, good restraint speckle is to the interference of magnetic resonance image (MRI), and it is rebuild image property and is improved and improves.

As preferred method, step 2 specifically comprises the following steps:

A. based on image area method, equation description is carried out to coil aliased image I: I=E ρ;

Wherein, E is coil sensitivity matrix, and ρ is image to be reconstructed;

B. utilize least squares estimate to solve ρ, namely minimize following objective function:

$\mathrm{\ρ}=\arg min\left\{\right||E\mathrm{\ρ}-I|{|}_2^2\}$

C. introduce regularization method, obtain following formula;

ρ _λ＝arg min{||Eρ-I|| ₂+λ ²||L(ρ-ρ ₀)|| ₂}

Wherein, λ is regularization parameter, and L is regular operator, ρ ₀for prior imformation;

D. utilize described L-curve method to formula ρ _λ=arg min{||E ρ-I|| ₂+ λ ²|| L (ρ-ρ ₀) || ₂calculate, obtain making model error || E ρ-I|| ₂with prior imformation error λ ²|| L (ρ-ρ ₀) || ₂sum reaches minimum optimum regularization parameter, obtains as shown in the formula ρ=(E ^he+ λ ²l ^hl) ^-1(E ^hi+ λ ²l ^hl ρ ₀);

E. by formula ρ=(E ^he+ λ ²l ^hl) ^-1(E ^hi+ λ ²l ^hl ρ ₀) be deformed into as shown in the formula (E ^he+ λ ²l ^hl) ρ=E ^hi+ λ ²l ^hl ρ ₀.

Adopting above-mentioned preferred scheme, preparing for calculating each coil restructured images.

As preferred method, described prior imformation ρ ₀it is the low-resolution image obtained by K space center data.

Adopt above-mentioned preferred scheme, introduce prior imformation, be convenient to calculate optimum regularization parameter.

As preferred method, step 3 is specially: utilize conjugate gradient iterative procedure to formula (E ^he+ λ ²l ^hl) ρ=E ^hi+ λ ²l ^hl ρ ₀carry out iterative processing, until error is less than the precision of setting, last iteration result is the reconstruction image of coil.

Adopt above-mentioned preferred scheme, effectively can obtain the image after rebuilding, resolution is high, and robustness is good.

Accompanying drawing explanation

Fig. 1 is the flow process frame diagram of a kind of parallel MR image rebuilding method based on regularization iteration of the present invention;

Fig. 2 is least square method, the reconstructed results of the inventive method and corresponding Error Graph (R=2);

A () is reference picture;

B reconstruction image that () obtains for adopting least square method;

(b1) the corresponding Error Graph for adopting least square method to obtain;

C reconstruction image that () obtains for adopting method of the present invention;

(c1) the corresponding Error Graph for adopting method of the present invention to obtain;

Fig. 3 is least square method, the reconstructed results of the inventive method and corresponding Error Graph (R=4);

A () is reference picture;

B reconstruction image that () obtains for adopting least square method;

(b1) the corresponding Error Graph for adopting least square method to obtain;

C reconstruction image that () obtains for adopting method of the present invention;

(c1) the corresponding Error Graph for adopting method of the present invention to obtain;

Fig. 4 for adding human brain data reconstruction result after spike noise in loop data;

A () is reference picture;

B reconstruction image that () obtains for adopting least square method;

C reconstruction image that () obtains for adopting method of the present invention.

Embodiment

The preferred embodiment of the present invention is described in detail below in conjunction with accompanying drawing.

In order to reach object of the present invention, in a kind of some of them embodiment of the parallel MR image rebuilding method based on regularization iteration,

As shown in Figure 1, based on the parallel MR image rebuilding method of regularization iteration, method comprises the following steps:

1. pair full K space data carries out analog sampling, obtains the K space data of the capable data of multi-channel parallel coil K space center self-calibration ACS and lack sampling, calculates each coil space sensitivity profile and each coil aliased image respectively;

2. in the computation process that coil aliased image is launched, in conjunction with least squares theory, based on image area method, equation description is carried out to reconstruction image, introduce regularization method, utilize L-curve method to calculate optimum regularization parameter, obtain the matrix equation after regularization;

3. utilize conjugate gradient iterative procedure to carry out iterative approximation to the matrix equation after regularization, obtain the reconstruction image of each coil.

The present invention proposes a kind of parallel MR image rebuilding method based on regularization iteration, the method is based on Least Square Theory, introduce regularization method in the process calculated, utilize L-curve method to solve optimum regularization parameter, carry out image reconstruction in conjunction with conjugate gradient iterative procedure.While improving signal to noise ratio (S/N ratio), strengthen the robustness of the method, effectively solve the problems of the technologies described above, reduce noise to large extent to the impact of reconstructed results, improve signal noise ratio (snr) of image.

In order to better reflect the superiority of the parallel MR image rebuilding method of a kind of regularization iteration of the present invention, first the parallel MR imaging method based on image area is described.

Suppose in magnetic resonance imaging process, receiving coil for having n receiving coil cellular array coil, the aliased image pixel value I that i-th receiving coil obtains _i(x, y) is as shown in the formula (1);

I _i(x，y)＝C _i(x，y ₁)ρ(x，y ₁)+C _i(x，y ₂)ρ(x，y ₂)+...+C _i(x，y _R)ρ(x，y _R) (1)

Wherein, C _ibe the plural sensitivity function of i-th coil unit, ρ is image to be reconstructed, and R accelerates multiple;

Formula (1) is converted, obtains as shown in the formula (2) and formula (3);

I＝Eρ (2)

$E=e^{{\mathrm{ik}}_{\mathrm{\κ}}{r}_p}{C}_j\left(r_p\right)---\left(3\right)$

Wherein, E is coil sensitivity matrix, r _pp image pixel, k _κfor a kth K spatial value, C _j(r _p) be positioned at r for a jth receiving coil _pthe coil sensitivity functions at place.

To formula (2), when known I, utilize the method for least-squares estimation to calculate ρ, namely minimize as shown in the formula (7):

$\mathrm{\ρ}=\arg min\left\{\right||E\mathrm{\ρ}-I|{|}_2^2\}---\left(7\right)$

For formula (7), the solution under signal to noise ratio (S/N ratio) optimal situation is formula (8):

ρ＝(E ^Hψ ^-1E) ^-1E ^Hψ ^-1I (8)

Wherein: E ^hbe the associate matrix of coil sensitivity matrix E, ψ is the noise correlation matrix of receiving coil unit, when ignoring noise, is usually taken as unit matrix.

By above formula (8), the operation of separating aliasing is carried out to each pixel in image area, just can reconstruct target image.

Formula (2) is an over-determined systems in fact, and the conditional number Cond (E) >=1 of its coil sensitivity matrix E, therefore the Solve problems of this system of linear equations is ill, and its solution has ill-posedness.But be often subject to the impact of motion and noise in the multi channel imaging data acquisition coil data of reality and produce exception, and destroy the situation that data and exceptional value often cause large residual error, and least square solution is usually very responsive to the data point with large residual error, cause separating instability, this will cause the bigger error of the least square estimation method, makes existing method not possess robustness.

If solely carry out pMRI image reconstruction by least square method, can be deteriorated to the robustness of the exceptional value such as artifact and noise.Be a selection to addressing this problem regularization, but the method is used for parallel MR imaging result is not also desirable, reason is that regularizing functionals is excessively smooth, causes the loss in detail of rebuild image.In addition, in all multi-methods of MR image reconstruction, introduce iterative approximation and not only reduce the computing time of computing machine when separating large-scale linear problem and requirement to internal memory, and undertaken rebuilding by process of iteration and can obtain the higher image of precision, and itself there is inherent regularization effect.But the regularization effect of this inherence is not sufficient to the pathosis of compensated linear problem.

By analysis above, the present invention adopts the method (Regularized reconstruction using Conjugate Gradient iterations, ReCG) of conjugate gradient iterative procedure compute matrix equation regularization.The method process is: least square method improved, Tikhonov regularization is introduced to formula (2), optimum regularization parameter is obtained with L-curve, obtain the matrix equation after regularization, again conjugate gradient iterative approximation is carried out to it, thus show that resolution is higher, the good reconstructed results of robustness.

Lucky big vast promise husband (Tikhonov) method of regularization is a kind of effective ways solving ill-posed problem, and its main thought introduces the prior imformation of separating to former problem, obtains and make the weighted array of the norm of solution and residual norm be minimum solution.

In the computation process of coil without aliased image, introduce regularization method, formula (2) is converted to as shown in the formula (4);

ρ _λ＝arg min{||Eρ-I|| ₂+λ ²||L(ρ-ρ ₀)|| ₂} (4)

Wherein, λ is regularization parameter, and L is regular operator, ρ ₀prior imformation, prior imformation ρ ₀it is the low-resolution image obtained by K space center data.

L regular operator is generally high-pass filtering operator.

Utilize L-curve method to calculate formula (4), obtain making model error || E ρ-I|| ₂with prior imformation error λ ²|| L (ρ-ρ ₀) || ₂sum reaches minimum optimum regularization parameter, by formula (4) through conversion, obtains following formula (5);

ρ＝(E ^HE+λ ²L ^HL) ^-1(E ^HI+λ ²L ^HLρ ₀) (5)

In order to avoid the calculating internal memory caused because of matrix inversion is occupied, formula (5) is converted to formula (6);

(E ^HE+λ ²L ^HL)ρ＝E ^HI+λ ²L ^HLρ ₀(6)

Next, utilize conjugate gradient iterative procedure to carry out iterative processing to formula (6), last iteration result is the reconstruction image of coil.Conjugate gradient iterative procedure to solution large linear systems and Large Scale Nonlinear optimization problem very effective, memory space is little, and convergence is high, good stability.

For the feasibility of a kind of parallel MR image rebuilding method based on regularization iteration of checking the present invention, the experimental data adopted is the full K space data (Gradient echoes sequence of 8 passage brain that the GE magnetic resonance scanning system of 1.5T obtains, TE/TR=4.1/9.7ms, T1=450ms, BW=15kHz, matrix is 200 × 200).Complete K space data generating reference image, as a comparison.The emulation of data lack sampling is carried out by the mode extracting some phase encode line, near retaining K space center, 24 row low-frequency datas are for estimating each coil space sensitivity profile, periphery is 2,4 carry out lack sampling with speedup factor R respectively, obtain the K space data lacked, carry out numerical value emulation method again and image reconstruction is carried out to method of the present invention, and carry out qualitative, quantitative contrast with the image reconstruction quality of least square method.

The MATLAB of simulation calculation on i5-4570CPU, 4GB memory computer carries out.

At the experiment initial stage, first K space data is emulated, draw reconstructed results by the method that least square method and the present invention propose respectively.When Fig. 2,3 represents speedup factor R=2, R=4 respectively, reference picture and two kinds of method reconstructed results and corresponding Error Graph.

The true brain data reconstructed results of Fig. 2 and Fig. 3 reaches Expected Results, by finding out with the qualitative contrast of original method, the method that the present invention proposes is to the reconstruction image that the suppression of noise and image detail resolution are all better than being obtained by least square method, when speedup factor is larger, effect is more obvious.Error image illustrates method reconstructed results reconstruction precision partially and method that the present invention proposes more intuitively to effective suppression of noise, more close to reference picture in intensity profile.

In order to embody method of the present invention further to the rejection ability of the exceptional values such as noise, then destruction data are emulated, add spike noise (density is 0.01) wherein in a loop data, give the reconstructed results of two kinds of methods during R=2, carry out contrast experiment.

Fig. 4 is the result of rebuilding loop data, and wherein (a) carries out the reference picture that sum-of-squares obtains for raw k-space data; (b) image for adopting the K space data destroyed to carry out least square method reconstruction.Can clearly see in the picture, in single parallel coil, the destruction of data can impact final reconstruction quality, and in reconstructed results, artifact is obvious.C reconstruction image that () obtains for adopting method of the present invention.Compare the image in (b), the interference that loop data noise causes almost is eliminated completely.In addition, what the details resolution of image was also expressed is very clear.

In addition, not single from visual angle, the image that method intensity profile adopting the present invention propose rebuilds out is also closer to reference picture.The present invention, when the quality of image is rebuild in quantitative test, adopts normalized mean squared error (NMSE) and signal to noise ratio (S/N ratio) (SNR) two indexs as image quality evaluation standard.Normalized mean squared error describes the relative error of rebuilding image and standard reference image, and computing method are as shown in the formula (9);

$NMSE=\frac{\underset{(x,y)}{\Σ}\left|\right|I_{recon}(x,y)|-|I_{ref}(x,y)|{|}^2}{\underset{(x,y)}{\Σ}{\left|I_{ref}(x,y)\right|}^2}---\left(9\right)$

Wherein: I _recon(x, y) image for reconstructing, I _ref(x, y) is reference picture.

Signal to noise ratio (S/N ratio) is also an importance of evaluation image quality, and its value is larger, and show that the noise spot in image is fewer, the quality of image is higher.The method of the study general employing of parallel MR imaging is as shown in the formula (10);

SNR＝10lg(S/N) (10)

S＝var(I _ref(：)，1) (11)

N＝mean(I _recon(：)-I _ref(：))^2(12)

S is the covariance of standard picture, illustrates the useful information of image, and N represents standard picture and rebuilds the square error of image, represents the noisiness of rebuilding image.

The signal to noise ratio (S/N ratio) calculated under different speedup factor condition and normalized mean squared error are as shown in Table 1 and Table 2.Can find out that from table 1 the inventive method obtains the higher reconstruction image of signal to noise ratio (S/N ratio), particularly when accelerating multiple and being larger, signal to noise ratio (S/N ratio) is greatly improved.Table 2 adopts normalized mean squared error (NMSE) as estimating, and from numerically can clearly see, the inventive method has better rejection ability to noise, and reconstruction precision is more close to reference picture, and when speedup factor is larger, advantage is more obvious.

Table 1 two kinds of methods rebuild the comparison of signal to noise ratio (S/N ratio) (SNR)

Table 2 two kinds of methods rebuild the comparison % of normalized mean squared error (NMSE)

On the basis analyzing parallel MR image reconstruction problem pathosis, consider the impact of noise on reconstructed results, propose a kind of new method of rebuilding for parallel imaging based on regularization conjugate gradient iteration.For pathosis and the signal noise ratio (snr) of image decline problem of parallel MR imaging process, regularization method and process of iteration are combined, effectively improves the robustness of method, rebuild image to the rejection ability of noise.Shown by the simulation experiment result simultaneously, when speedup factor is larger, utilize method for reconstructing used in the present invention still can to large extent stress release treatment produce adverse effect, make rebuild image signal to noise ratio (S/N ratio) be improved significantly.

Above-described is only the preferred embodiment of the present invention, it should be pointed out that for the person of ordinary skill of the art, and without departing from the concept of the premise of the invention, can also make some distortion and improvement, these all belong to protection scope of the present invention.

Claims (4)

1., based on the parallel MR image rebuilding method of regularization iteration, it is characterized in that, said method comprising the steps of:

1. pair full K space data carries out analog sampling, obtains the K space data of the capable data of multi-channel parallel coil K space center self-calibration ACS and lack sampling, calculates each coil space sensitivity profile and each coil aliased image respectively;

2. in the computation process that coil aliased image is launched, in conjunction with least squares theory, based on image area method, equation description is carried out to reconstruction image, introduce regularization method, utilize L-curve method to calculate optimum regularization parameter, obtain the matrix equation after regularization;

3. utilize conjugate gradient iterative procedure to carry out iterative approximation to the matrix equation after regularization, obtain the reconstruction image of each coil.

2. the parallel MR image rebuilding method based on regularization iteration according to claim 1, it is characterized in that, described step 2 specifically comprises the following steps:

A. based on image area method, equation description is carried out to described coil aliased image I: I=E ρ;

Wherein, E is coil sensitivity matrix, and ρ is image to be reconstructed;

B. utilize least squares estimate to solve ρ, namely minimize following objective function:

$\mathrm{\ρ}=\arg min\left\{\right||E\mathrm{\ρ}-I|{|}_2^2\}$

C. introduce regularization method, obtain following formula;

${\mathrm{\ρ}}_{\mathrm{\λ}}=\arg min\left\{\right||E\mathrm{\ρ}-I|{|}_2+{\mathrm{\λ}}^2\left|\right|L(\mathrm{\ρ}-{\mathrm{\ρ}}_0)\left|{|}_2\right\}$

Wherein, λ is regularization parameter, and L is regular operator, ρ ₀for prior imformation;

D. utilize described L-curve method to formula ρ _λ=arg min{||E ρ-I|| ₂+ λ ²|| L (ρ-ρ ₀) || ₂calculate, obtain making model error || E ρ-I|| ₂with prior imformation error λ ²|| L (ρ-ρ ₀) || ₂sum reaches minimum optimum regularization parameter, obtains as shown in the formula ρ=(E ^he+ λ ²l ^hl) ^-1(E ^hi+ λ ²l ^hl ρ ₀);

E. by formula ρ=(E ^he+ λ ²l ^hl) ^-1(E ^hi+ λ ²l ^hl ρ ₀) be deformed into as shown in the formula (E ^he+ λ ²l ^hl) ρ=E ^hi+ λ ²l ^hl ρ ₀.

3. the parallel MR image rebuilding method based on regularization iteration according to claim 2, is characterized in that, described prior imformation ρ ₀it is the low-resolution image obtained by K space center data.

4. the parallel MR image rebuilding method based on regularization iteration according to claim 3, it is characterized in that, described step 3 is specially: utilize described conjugate gradient iterative procedure to formula (E ^he+ λ ²l ^hl) ρ=E ^hi+ λ ²l ^hl ρ ₀carry out iterative processing, until error is less than the precision of setting, last iteration result is the reconstruction image of coil.