CN106108903A - 一种改进的并行磁共振图像重建方法 - Google Patents

Abstract

本发明公开了一种改进的并行磁共振图像重建方法，属于磁共振成像领域，适用于并行磁共振成像过程中K空间数据的重建方法。本方法在拟合权重系数时，利用欠采样数据附近的K空间数据区域与自校准数据区域进行估算，而不是利用邻近点和自校准线，如此可以较好的拟合出K空间数据的非线性关系；在奇异值分解并进行截尾处理时可以抑制采样信号中的噪声，保证了图像信噪比较高。本发明的方法是基于矩阵广义逆和奇异值分解的GRAPPA改进方法，可抑制K空间数据的噪声，提高拟合精确度，同时在加速因子较大时仍能保证成像质量。

Description

一种改进的并行磁共振图像重建方法

技术领域

本发明属于磁共振成像领域，涉及并行磁共振成像方法，具体地说是一种适用于并行磁共振成像过程中K空间数据的重建方法。

背景技术

近几十年，磁共振成像(Magnetic Resonance Imaging，MRI)一直是医学成像的热点问题，被广泛用于临床诊断和科学研究等领域。传统的磁共振成像对K空间数据进行全采样，成像速度较慢，且成像质量容易受到呼吸、血液流动等运动因素影响而出现伪影。针对传统磁共振成像中出现的问题，人们提出了许多并行磁共振成像技术。例如在文献：Griswold MA,Jakob PM,Heidemann RM,et al.Generalized autocalibrating partiallyparallel acquisition.Magn Reson Med.2002；47:1202-1210中，提出的GRAPPA(Generalized autocalibrating partially parallel acquisition,GRAPPA)算法，对K空间数据进行欠采样，利用ACS进行拟合得到权重系数，将权重系数应用到K空间的其他欠采样区域，由此得到完整的K空间数据。但是，GRAPPA算法拟合的准确度过于依赖权重系数，实际采集数据中的噪声将对重建结果产生较大的偏差。

发明内容

针对现有的GRAPPA算法的不足，本发明提出一种基于矩阵广义逆和奇异值分解的GRAPPA改进方法。可抑制K空间数据的噪声，提高拟合精确度，同时在加速因子较大时仍能保证成像质量。

用于实现本发明的技术解决方案包括如下步骤：

一种改进的并行磁共振图像重建方法，包括如下步骤：

步骤1：对成像物体进行欠采样，得到数据矩阵D，全采样校准区域，得到数据阵F；

步骤2：对数据矩阵D进行数据提取和重新构造，得到数据矩阵B；

步骤3：对数据矩阵F进行数据提取和重新构造，得到数据矩阵A；

步骤4：对数据矩阵A进行奇异值分解；

步骤5：对步骤4中得到的奇异值进行截尾处理；

步骤6：设置计数变量t＝1；

步骤7：检测数据矩阵B的第t行数据中的非零元素，将所有非零元素下标记为集合E_t，非零元素个数记为e_t；变量下面加下标t表示的是与第t次递归有关的变量；

步骤8：利用截尾处理的数据矩阵A_t，抽值得到矩阵

步骤9：计算矩阵的Moore-Penrose逆矩阵

步骤10：重建出数据矩阵B的第t行数据r_t，并使用向量r_t替代B(t,:)；

步骤11：判断计数变量t是否等于S_x×S_y，如果不等于，则将计数变量t加1，并返回步骤7，否则进入步骤12；

步骤12：重新设置计数变量q＝1；

步骤13：利用数据矩阵B的列向量B(:,q)重新构造成S_x×S_y维数据矩阵

步骤14：将矩阵扩展成P×Q×C维的数据矩阵其中数据区域 的其他区域用0进行填充；

步骤15：判断计数变量q是否等于L²，如果不等于，则将计数q加1，并返回步骤13，否则进入步骤16；

步骤16：重建出K空间数据矩阵其中i＝1,2,…L，j＝1,2,…,L，h＝1，2，…C，用D中数值不为0的数据替代R中相应位置的数据；

步骤17：由重建数据矩阵R进行傅里叶逆变换得到各线圈的K空间数据，对数据进行傅里叶逆变换得到各线圈的图像，并对各线圈图像求平方和，得到最终的重建结果。

进一步地，所述步骤1的具体过程包括：

相位编码方向采用欠采样扫描方式，频率编码方向采用全采样扫描方式，获得K空间数据记为P×Q×C维的数据矩阵D，其中P表示相位编码方向的数据线数目，Q表示频率编码方向的数据线数目，C表示扫描线圈的个数；K空间数据的中心区域为校准区域，其采用全采样方式扫描，若校准区域在相位编码方向的维数为M，则该区域对应的数据矩阵其中P、M和Q的取值均设置为偶数。

进一步地，所述步骤2的具体方法包括：

使用S_x×S_y维的扫描窗口对数据矩阵D进行数据提取和重新构造，扫描次数为N，其中S_x＝P-L+1，S_y＝Q-L+1，N＝L²C，L的取值人为设定且满足(M-L+1)S_y≥N，数据提取、重新构造的步骤如下：

步骤2.1：以d_i,j,h为起点对矩阵D进行数据扫描，得到数据矩阵W_i,j,h，其中d_i,j,h＝[D]_i,j,h，i＝1,2,…,L，j＝1,2,…,L，h＝1,2,…,C；

步骤2.2： 为数据矩阵W_i,j,h的列向量：

$w_{i,j,h}^p=\left[\begin{array}{c}{d}_{i,j+p-1,h}\\ d_{i+1,j+p-1,h}\\ .\\ .\\ .\\ d_{i+S_x-1,j+p-1,h}\end{array}\right]$

其中p＝1,2,…,S_y；

步骤2.3：利用构造S_x×S_y维列向量m_i,j,h：

$m_{i,j,h}=\left[\begin{array}{c}{w}_{i,j,h}^1\\ w_{i,j,h}^2\\ .\\ .\\ .\\ w_{i,j,h}^{S_y}\end{array}\right];$

步骤2.4：利用m_i,j,h构造出S_xS_y×N维数据矩阵其中 i_g＝g-L²×(h_g-1)-L×(j_g-1)

表示向下取整，g＝1,2,…,N。

进一步地，所述步骤3的具体实现方法与所述步骤2的具体实现方法相同。

进一步地，所述步骤4中对数据矩阵A进行奇异值分解的具体实现包括：

$A=\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\σ}}_k{u}_k{v}_k^H$

其中(·)^H表示共轭转置，u_k为A的左奇异值向量，v_k为A的右奇异值向量，σ_k为矩阵A的奇异值。

进一步地，所述步骤8的具体过程包括：

步骤8.1：构造矩阵z_t＝min(z,e_t)，z为截尾处理后不为零的奇异值的个数；定义为A_t的第n行向量；

步骤8.2：构造校准矩阵 其中

进一步地，所述步骤13的数据矩阵具体为：

$R_{i_r,j_r,h_r}=\left[\begin{array}{cccc}{b}_{1,q}& b_{S_x+1,q}& \mathrm{...}& b_{(S_y-1)S_x,q}\\ b_{2,q}& b_{S_x+2,q}& \mathrm{...}& b_{(S_y-1)S_x+1,q}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ b_{S_x,q}& b_{2S_x,q}& \mathrm{...}& b_{S_y{S}_x,q}\end{array}\right]$

其中i_r＝q-L²×(h_r-1)-L×(j_r-1)，i_b＝1,2,…,S_xS_y，j_b＝1，2，…，N。

本发明的有益效果：

1、本方法在拟合权重系数时，利用欠采样数据附近的K空间数据区域与自校准数据区域进行估算，而不是利用邻近点和自校准线。因此，可以较好的拟合出K空间数据的非线性关系。

2、本方法在奇异值分解并进行截尾处理时可以抑制采样信号中的噪声，保证了图像信噪比较高。

附图说明

图1是本发明的实施流程图。

图2是仿真实验的全K空间参考图像。

图3是本发明的重建图像。

图4是GRAPPA算法的重建图像。

图5是本发明的重建图像与参考图像的幅度差值图像。

图6是GRAPPA算法的重建图像与参考图像的幅度差值图像。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

如图1所示，本发明实现的具体流程图如下：

(1)相位编码方向采用欠采样扫描方式，频率编码方向采用全采样扫描方式，获得K空间数据记为P×Q×C维的数据矩阵D，其中P表示相位编码方向的数据线数目，Q表示频率编码方向的数据线数目，C表示扫描线圈的个数。K空间数据的中心区域为校准区域，其采用全采样方式扫描。若校准区域在相位编码方向的维数为M，则该区域对应的数据矩阵其中P、M和Q的取值均假定为偶数。

(2)使用S_x×S_y维的扫描窗口对数据矩阵D进行数据提取和重新构造，扫描次数为N，其中S_x＝P-L+1，S_y＝Q-L+1，N＝L²C，L的取值人为设定且满足(M-L+1)S_y≥N，数据提取、重新构造的步骤如下：

以d_i,j,h为起点对矩阵D进行数据扫描，得到数据矩阵W_i,j,h，其中d_i,j,h＝[D]_i,j,h，i＝1,2,…,L，j＝1,2,…,L，h＝1,2,…,C。

为数据矩阵W_i,j,h的列向量：

$w_{i,j,h}^p=\left[\begin{array}{c}{d}_{i,j+p-1,h}\\ d_{i+1,j+p-1,h}\\ .\\ .\\ .\\ d_{i+S_x-1,j+p-1,h}\end{array}\right]$

其中p＝1,2,…,S_y。

利用构造S_x×S_y维列向量m_i,j,h：

$m_{i,j,h}=\left[\begin{array}{c}{w}_{i,j,h}^1\\ w_{i,j,h}^2\\ .\\ .\\ .\\ w_{i,j,h}^{S_y}\end{array}\right].$

利用m_i,j,h构造出S_xS_y×N维数据矩阵其中

i_g＝g-L²×(h_g-1)-L×(j_g-1)

表示向下取整，g＝1,2,…,N。

(3)使用S'_x×S_y维的扫描窗口对数据矩阵F进行数据扫描和重新构造得到S'_xS_y×N维数据矩阵A，扫描次数为N，其中S'_x＝M-L+1，S_y＝Q-L+1，N＝L²C，数据提取和构造的方式同步骤(2)类似。

(4)对A进行奇异值分解：

$A=\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\σ}}_k{u}_k{v}_k^H$

其中(·)^H表示共轭转置，u_k为A的左奇异值向量，v_k为A的右奇异值向量，σ_k为矩阵A的奇异值。

(5)对矩阵A的奇异值进行截尾处理：若则σ_k＝0，k＝1,2,...,N，本专利中取ε＝0.05，记z为截尾处理后不为零的奇异值的个数。

(6)设置计数变量t＝1。

(7)检测数据矩阵B的第t行数据中的非零元素，将所有非零元素下标记为集合E_t，非零元素个数记为e_t；变量下面加下标t表示的是与第t次递归有关的变量。

(8)构造矩阵z_t＝min(z,e_t)。定义为A_t的第n行向量，进一步构造校准矩阵 其中

(9)计算矩阵的Moore-Penrose逆矩阵

(10)利用矩阵B的第t行数据B(t,:)重建出N维行向量r_t：使用向量r_t替代B(t,:)。

(11)判断计数变量t是否等于S_x×S_y，如果不等于，则将计数变量t加1，并返回步骤(7)，否则重建结束，进入步骤(12)。

(12)重新设置计数变量q＝1。

(13)利用数据矩阵B的列向量B(:,q)重新构造成S_x×S_y维数据矩阵

$R_{i_r,j_r,h_r}=\left[\begin{array}{cccc}{b}_{1,q}& b_{S_x+1,q}& \mathrm{...}& b_{(S_y-1)S_x,q}\\ b_{2,q}& b_{S_x+2,q}& \mathrm{...}& b_{(S_y-1)S_x+1,q}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ b_{S_x,q}& b_{2S_x,q}& \mathrm{...}& b_{S_y{S}_x,q}\end{array}\right]$

其中i_b＝1,2,…,S_xS_y，j_b＝1，2，…，N, i_r＝q-L²×(h_r-1)-L×(j_r-1)。

(14)将矩阵扩展成P×Q×C维数据矩阵其中数据区域 的其他区域用0进行填充。

(15)判断计数变量q是否等于L²，如果不等于，则将计数q加1，并返回步骤(13)，否则重构结束，进入步骤(16)。

(16)重建数据矩阵其中i＝1,2,…L，j＝1,2,…,L，h＝1，2，…C。

用D中数值不为0的数据替代R中相应位置的数据。

(17)由重建数据矩阵R可以得到各线圈的K空间数据，对数据进行傅里叶逆变换得到各线圈的图像，并对各线圈图像求平方和，得到最终的重建结果。

下面结合仿真实验对本发明的效果做进一步说明。

为了评估本方法的性能，考虑一人体大脑K空间数据的重建过程，加速因子R＝5，线圈个数为8，背景噪声假设为高斯白噪声。

实验条件

实验1，分别使用本发明和GRAPPA算法对K空间数据进行重建，重建后图像分别为图3和图4，图2为全K空间参考图像。

实验2，将本发明和GRAPPA算法得到的重建图像的幅度值分别与全K空间参考图像幅度值相减，得到幅度值误差图像分别如图5和图6所示。

实验分析

从图2、图3、图4的相互对比可以看出本算法与GRAPPA算法相比，噪声和伪影明显较少，图像纯净度较高，较为接近参考图像。

从图5，图6的对比可以看出本发明重建图像较为接近全K空间参考图像，本发明重建图像与全K空间参考图像方差值明显小于GRAPPA算法与全K空间参考图像方差值。本方法重建图像较为接近参考图像。

上文所列出一系列的详细说明并非用以限制本发明的保护范围，凡未脱离本发明技艺精神所作的等效实施方式或变更均应包含在本发明的保护范围之内。

Claims (7)

1.一种改进的并行磁共振图像重建方法，其特征在于，包括如下步骤：

步骤1：对成像物体进行欠采样，得到数据矩阵D，全采样校准区域，得到数据阵F；

步骤2：对数据矩阵D进行数据提取和重新构造，得到数据矩阵B；

步骤3：对数据矩阵F进行数据提取和重新构造，得到数据矩阵A；

步骤4：对数据矩阵A进行奇异值分解；

步骤5：对步骤4中得到的奇异值进行截尾处理；

步骤6：设置计数变量t＝1；

步骤7：检测数据矩阵B的第t行数据中的非零元素，将所有非零元素下标记为集合E_t，非零元素个数记为e_t；变量下面加下标t表示的是与第t次递归有关的变量；

步骤8：利用截尾处理的数据矩阵A_t，抽值得到矩阵

步骤9：计算矩阵的Moore-Penrose逆矩阵

步骤10：重建出数据矩阵B的第t行数据r_t，并使用向量r_t替代B(t,:)；

步骤11：判断计数变量t是否等于S_x×S_y，如果不等于，则将计数变量t加1，并返回步骤7，否则进入步骤12；

步骤12：重新设置计数变量q＝1；

步骤13：利用数据矩阵B的列向量B(:,q)重新构造成S_x×S_y维数据矩阵

步骤14：将矩阵扩展成P×Q×C维的数据矩阵其中数据区域 的其他区域用0进行填充；

步骤15：判断计数变量q是否等于L²，如果不等于，则将计数q加1，并返回步骤13，否则进入步骤16；

步骤16：重建出K空间数据矩阵其中i＝1,2,…L，j＝1,2,…,L，

h＝1，2，…C，用D中数值不为0的数据替代R中相应位置的数据；

步骤17：由重建数据矩阵R进行傅里叶逆变换得到各线圈的K空间数据，对数据进行傅里叶逆变换得到各线圈的图像，并对各线圈图像求平方和，得到最终的重建结果。

2.根据权利要求1所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤1的具体过程包括：

相位编码方向采用欠采样扫描方式，频率编码方向采用全采样扫描方式，获得K空间数据记为P×Q×C维的数据矩阵D，其中P表示相位编码方向的数据线数目，Q表示频率编码方向的数据线数目，C表示扫描线圈的个数；K空间数据的中心区域为校准区域，其采用全采样方式扫描，若校准区域在相位编码方向的维数为M，则该区域对应的数据矩阵其中P、M和Q的取值均设置为偶数。

3.根据权利要求1所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤2的具体方法包括：

使用S_x×S_y维的扫描窗口对数据矩阵D进行数据提取和重新构造，扫描次数为N，其中S_x＝P-L+1，S_y＝Q-L+1，N＝L²C，L的取值人为设定且满足(M-L+1)S_y≥N，数据提取、重新构造的步骤如下：

步骤2.1：以d_i,j,h为起点对矩阵D进行数据扫描，得到数据矩阵W_i,j,h，其中d_i,j,h＝[D]_i,j,h，i＝1,2,…,L，j＝1,2,…,L，h＝1,2,…,C；

步骤2.2： 为数据矩阵W_i,j,h的列向量：

$w_{i,j,h}^p=\left[\begin{array}{c}{d}_{i,j+p-1,h}\\ d_{i+1,j+p-1,h}\\ \begin{array}{c}.\\ .\\ .\end{array}\\ d_{i+S_x-1,j+p-1,h}\end{array}\right]$

其中p＝1,2,…,S_y；

步骤2.3：利用构造S_x×S_y维列向量m_i,j,h：

$m_{i,j,h}=\left[\begin{array}{c}{w}_{i,j,h}^1\\ w_{i,j,h}^2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ w_{i,j,h}^{S_y}\end{array}\right];$

步骤2.4：利用m_i,j,h构造出S_xS_y×N维数据矩阵其中 i_g＝g-L²×(h_g-1)-L×(j_g-1)

表示向下取整，g＝1,2,…,N。

4.根据权利要求3所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤3的具体实现方法与所述步骤2的具体实现方法相同。

5.根据权利要求1所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤4中对数据矩阵A进行奇异值分解的具体实现包括：

$A=\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\σ}}_k{u}_k{v}_k^H$

其中(·)^H表示共轭转置，u_k为A的左奇异值向量，v_k为A的右奇异值向量，σ_k为矩阵A的奇异值。

6.根据权利要求1所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤8的具体过程包括：

步骤8.1：构造矩阵z_t＝min(z,e_t)，z为截尾处理后不为零的奇异值的个数；定义为A_t的第n行向量；

步骤8.2：构造校准矩阵 其中

7.根据权利要求1所述的一种改进的并行磁共振图像重建方法，其特征在于，所述步骤13的数据矩阵具体为：

$R_{i_r,j_r,h_r}=\left[\begin{array}{cccc}{b}_{1,q}& b_{S_x+1,q}& ...& b_{(S_y-1)S_x,q}\\ b_{2,q}& b_{S_x+2,q}& ...& b_{(S_y-1)S_x+1,q}\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ b_{S_x,q}& b_{2S_x,q}& ...& b_{S_y{S}_x,q}\end{array}\right]$

其中i_r＝q-L²×(h_r-1)-L×(j_r-1)，i_b＝1,2,…,S_xS_y，j_b＝1，2，…，N。